function x = forbacksubs(A,x,f)

% FORWARD (f=0). Forward elimination.
% For lower triangular L, x = forbacksubs(L,b,0) solves L*x = b.
% BACKSUBS (f=1).  Back substitution.
% For upper triangular U, x = forbacksubs(U,b,1) solves U*x = b.


[n,n] = size(A);
if f==0
    x(1) = x(1)/A(1,1);
    for k = 2:n
        j = 1:k-1;
        x(k) = (x(k) - A(k,j)*x(j))/A(k,k);
    end
    
else
    
    x(n) = x(n)/A(n,n);
    
    for k = n-1:-1:1
        j = k+1:n;
        x(k) = (x(k) - A(k,j)*x(j))/A(k,k);
    end
end
