function [L,U] =lupiv(A)
%LU-factorization with pivoting
[n,m]=size(A);
b=1:n;
I=eye(n);
LL=I;
P=ones(5,1)*[1:n];
if n~=m, disp('Your matrix is not square')
else
    for j=1:n-1
        %k=find(A(j:end,j)==max(A(j:end,j)));
        %k=k(1);
        %P(j,j)=k;P(j,k)=j;
        %p=P(j,:);
        %Perm=I(:,p);
        %B=Perm*A;
        B=A;
        L=I;
        L(j+1:end,j)=-B(j+1:end,j)/B(j,j);
        A=L*B;
        LL=LL/L;
    end
    L=LL;
    U=A;
end


end

