A=[4,1,-1,0;1,3,-1,0;-1,-1,5,2;0,0,2,4];
b=[7,8,-4,6]'
[m,n]=size(A);
l=eye(m);
r=zeros(m,n);
r(1,:)=A(1,:);
l(:,1)=A(:,1)/r(1,1);
for i=2:n
for j=i:n
r(i,j)=A(i,j)-l(i,1:i-1)*r(1:i-1,j);
end
for j=i+1:n
l(j,i)=(A(j,i)-l(j,1:i-1)*r(1:i-1,i))/r(i,i);
end
end
l
r
%l*y=b,A*x=b=>l*r*x=b
y(1)=b(1)/l(1,1);
for i=2:n
y(i)=(b(i)-l(i,1:i-1)*y(1:i-1)')/l(i,i)
end
%r*x=y
x(n)=y(n)/r(n:n);
for i=n-1:-1:1
x(i)=(y(i)-r(i,i+1:n)*x(i+1:n)')/r(i:i);
end
x