function y=lsolve(L,b);
%
% Given Ax=b and A = LU, solve Ly = b and return y
n = length(b); % length of the vector b
for i = 1:n
y(i) = b(i);
for j = 1:i-1 % solve for each y(i) using previously computed y(j), j = 1,...,i-1
y(i) = y(i) - L(i,j)*y(j);
end
end
y = y(:);