function [L,U]=lufactor(A)
% LU factorization for n x n matrix A
% L is lower triangular
% U is upper triangular
% A= L*U
% WARNING: possible division by zero

n=length(A);
U=A;
L=eye(n);
for p=1:n
     for q=p+1:n
           L(q,p)=U(q,p)/U(p,p);
           U(q,p:n)=U(q,p:n)-U(p,p:n)*L(q,p);
     end
end
