% Solving Ax=b by LU Factorization without pivoting
function [x,L,U]=LUsolve(A,B)
NA=size(A,2);
[NB1,NB]=size(B);
if NB1~=NA
    error('A and B must have compatible dimensions')
end
AB=[A B];
E=eye(NA);
L=eye(NA);
for j=1:NA-1
    for i=j+1:NA
        if (AB(j,j)~=0)
            E(i,j)=-AB(i,j)/AB(j,j);
            AB=E*AB;
            L=L*inv(E);
            E=eye(NA);
        else 
            error ('algorithm needs pivoting')
        end 
    end
end
U=AB(:,1:NA);
x=(U\(L\B));
