 function [L,U] = lufactor (A) 
%---------------------------------------------
%  function [L,U] = lufactor (A) 
%  performs LU factorization of a matrix A 
%---------------------------------------------
 n = size(A,1) ;
 for k=1:n-1 
     for i=k+1:n
       piv =  A(i,k) / A(k,k) ;
       A(i,k) = piv; 
       A(i,k+1:n) = A(i,k+1:n) - piv*A(k,k+1:n);
   end
   A
   pause
 end 
 L = eye(n) + tril(A,-1);
 U = triu(A);
 