%%
clear
n = 5;
format compact
%%
% It's easy to get just the lower triangular part of any matrix using the tril command.
A = magic(n)
L = tril(A)

%%
% We'll set up and solve a linear system with this matrix.
b = ones(n, 1);
x = forwardsub(L, b)

%%
% It's not clear what the error in this answer is. However, the
% residual, while not zero, is virtually $\epsilon_M$ in size.
b - L*x

%%
% It's easy to get just the upper triangular part of any matrix using the triu command.
A
U = triu(A)

%%
% We'll set up and solve a linear system with this matrix.
y = backsub(U, b)

%%
% It's not clear what the error in this answer is. However, the
% residual, while not zero, is virtually $\epsilon_M$ in size.
b - U*y