function [x] = lsol(L,b);
%
% solve the system Lx=b, where L is assumed to be unit lower
% triangular (the code looks only at the strictly lower triangular
% elements).