--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/control/lyap.m	Tue Aug 10 22:15:36 1993 +0000
@@ -0,0 +1,78 @@
+function x = lyap(a,b,c)
+
+% usage: x = lyap(a,b{,c})
+%
+% If (a,b,c) are specified, then lyap returns the solution of the Sylvester 
+%   equation
+%	a x + x b + c = 0
+%
+% If only (a,b) are specified, then lyap returns the solution of the 
+%   Lyapunov equation
+% 	a' x + x a + b = 0
+%   If b is not square, then lyap returns the solution of either
+%	a' x + x a + b' b = 0 	
+%   or
+%    	a x + x a' + b b' = 0
+%   whichever is appropriate.
+
+if((nargin ~= 3) && (nargin ~= 2))
+  error("lyap: illegal number of arguments")
+  return
+endif
+
+if( (n = is_square(a)) == 0 )
+  error("lyap: a is not square");
+  return
+endif
+
+if(nargin == 2)
+disp('LYAP')
+a=a
+b=b
+  % transform Lyapunov equation to Sylvester equation form
+  if ( (m=is_square(b)) == 0)
+    if( (m=rows(b)) == n)
+disp('TYPE CONTROL')
+      % solve a x + x a' + b b' = 0
+      b = b*b';
+      a = a';
+    else 
+disp('TYPE OBSERVE')
+      % (try to ) solve a'x + x a + b' b = 0
+      m = columns(b);
+      b = b'*b;
+    endif
+    if(m ~= n)
+      error("lyap: a, b not conformably dimensioned");
+      return;
+    endif
+  endif
+
+  % set up sylvester equation
+  c = b;
+  b = a;
+  a = b'
+else 
+  % check dimensions
+  if( (m = is_square(b)) == 0)
+    error("lyap: b must be square in a sylvester equation");
+    return
+  endif
+  [n1,m1] = size(c);
+  if((n ~= n1) || (m ~= m1))
+    error("lyap: a,b,c not conformably dimensioned");
+  endif
+endif
+a
+b
+c
+
+% call octave built-in function
+x = syl(a,b,c);
+
+disp('LYAP ERR')
+a*x + x*b + c
+a=a
+b=b
+c=c
+endfunction