function [L,U] = mylu(A);

% LU without pivoting

n = size(A,1);  if (size(A,2) ~= n), error('A must be square.');  end;

L = eye(n); U = zeros(n,n);

for i=1:(n-1)
   if (abs(A(i,i)) < 1e-10), 
      error('Zero pivot.  A cannot be factored as LU.'); 
   end;
   U(i,i:n) = A(i,i:n);
   L((i+1):n,i) = A((i+1):n,i)/A(i,i); 
   A((i+1):n,(i+1):n) = A((i+1):n, (i+1):n) - L((i+1):n,i)*U(i,(i+1):n);
end;
if (abs(A(n,n)) < 1e-10), 
   error('Zero pivot.  A cannot be factored as LU.'); 
end;
U(n,n) = A(n,n);