function [L,U]=mylu(A)
%
% [L,U] = mylu(A)
%
% This function implements Gaussian elimination
% with no pivoting. The routine is not robust
% and can fail even for nonsingular A
n = size(A,1);
for k = 1:n-1
for i = k+1:n
A(i,k) = A(i,k)/A(k,k); % Get multipliers
for j = k+1:n
A(i,j) = A(i,j) - A(i,k)*A(k,j); % Perform row sums
end
end
end
% Extract lower and upper triangular parts into L and U matrices
% NOTE: You would not normally do this because of storage
% considerations; this is only for demonstration purposes.
U = triu(A);
L = tril(A,-1) + eye(n);