function x=mylu(A,b)
%Input  -A is an N*N matrix
%       -b is N*1 matrix
%Output -x is an N*1 matrix containing the solution to Ax=b
[N,N]=size(A);x=zeros(N,1);y=zeros(N,1);
U0=A;L0=eye(N);L0t=zeros(N);P0=eye(N);R=zeros(N,N);
for p=1:N-1
    %Find the pivot row for colum p
    [max1,j]=max(abs(U0(p:N,p)));
    %Interchange row p and j
    R(p,:)=1:N;P1=eye(N);
    c=P1(p,:);P1(p,:)=P1(j+p-1,:);P1(j+p-1,:)=c;
    d=R(p,p);R(p,p)=R(p,j+p-1);R(p,j+p-1)=d;
    U0=P1*U0;P0=P1*P0;
    if A(p,p)==0 
        'A is singular.No unique solution' 
        break
    end
    L1=zeros(N);
    for k=p+1:N  
        mult=U0(k,p)/U0(p,p);U0(k,p:N)=U0(k,p:N)-mult*U0(p,p:N);L1(k,p)=-mult;
    end
    L0t=L0t+L1;
end
L0t=L0t+eye(N);R(N,:)=1:N;
for p=1:N-1
    L1=eye(N);L1(:,p)=L0t(:,p);
    for k=p:N-1 P2=eye(N);
        for n=1:N P1(:,n)=P2(:,R(k+1,n));end 
        L1=P1*L1*(eye(N)/P1);
    end
    L0=L0*(eye(N)/L1);
end
%Solve for y
b=P0*b;y(1)=b(1);for k=2:N y(k)=b(k)-L0(k,1:k-1)*y(1:k-1);end
%Solve for x
x(N)=y(N)/U0(N,N);for k=N-1:-1:1 x(k)=(y(k)-U0(k,k+1:N)*x(k+1:N))/U0(k,k);end