--- a/src/DLD-FUNCTIONS/qz.cc
+++ b/src/DLD-FUNCTIONS/qz.cc
@@ -287,6 +287,9 @@
     return (fabs (p) >= 1 ? 1 : -1);
 }
 
+
+//FIXME: Matlab does not produce lambda as the first output argument.
+//       Compatibility problem?
 DEFUN_DLD (qz, args, nargout,
   "-*- texinfo -*-\n\
 @deftypefn  {Loadable Function} {@var{lambda} =} qz (@var{A}, @var{B})\n\
@@ -1236,3 +1239,28 @@
 
   return retval;
 }
+
+/*
+%!shared a, b, c
+%!  a = [1 2; 0 3];
+%!  b = [1 0; 0 0];
+%!  c = [0 1; 0 0];
+%!assert(qz (a,b), 1);
+%!assert(isempty (qz (a,c)));
+
+%% Exaple 7.7.3 in Golub & Van Loan
+%!test
+%! a = [ 10  1  2;
+%!        1  2 -1;
+%!        1  1  2];
+%! b = reshape(1:9,3,3);
+%! [aa, bb, q, z, v, w, lambda] = qz (a, b);
+%! sz = length(lambda);
+%! assert ( (a*v) (:, 1:sz), (b * v * diag ([lambda;0])) (:, 1:sz), 1e-14);
+%! assert ( (w'*a) (1:sz, :) , (diag ([lambda;0]) * w' * b) (1:sz, :), 1e-14);
+%! assert (q * a * z, aa, 1e-14);
+%! assert (q * b * z, bb, 1e-14);
+
+%% FIXME: Still need a test for third form of calling qz
+
+*/