function x = fwsolve(A, b)
% Solves Ax = b for x, assuming that A is a nonsingular lower 
% triangular matrix.

n = length(b);
x = b;

% add tests for A(i,i) close to 0!
x(1) = b(1) / A(1,1);
for i = 2:n
    x(i) = (b(i) - A(i,1:i-1)*x(1:i-1)) / A(i,i);
end

return