function x=MyLR(A,y)
tic
% A=LR
[L R y]=LRDecomp(A,y);

% Ax=y <=> LRx=y <=> Lb=y AND Rx=b
b = LSolve(L,y);
x = RSolve(R,b);
toc
end

function x=RSolve(R,b)

n = length(b);

for k=n:-1:1
    x(k) = b(k);
    for i=k+1:n
        x(k) = x(k) - x(i)*R(k,i);
    end
    x(k) = x(k)/R(k,k);
end

end

function b=LSolve(L,y)

n = length(y);

for k=1:n
    b(k)=y(k);
    for i=1:k-1
        b(k) = b(k) - L(k,i)*b(i);
    end
end

end



function [L R y]=LRDecomp(A,y)


n = length(y);
x = zeros(1,n);

R = A;
L = eye(n);
for k=1:n-1
    for i=k+1:n
        L(i,k) = R(i,k)/R(k,k);
        R(i,k) =0;
        for j=k+1:n
            R(i,j) = R(i,j) - L(i,k)*R(k,j);
        end
    end
end



end

