function y = fs(L,b)
% FS  forward substitution with no fancy stuff
% usage: to solve the system  Ly = b,  do
%   >> x = fs(L,b) 

n = size(L,1);
y = zeros(n,1);

for j = 1:n
  sum = 0;
  for k = 1:j-1
    sum = sum + L(j,k) * y(k);
  end
  y(j) = (b(j) - sum) / L(j,j);
end
