function x = luab(A, b)
% Function to solve Ax = b using LU decomposition
%
% A is nxn matrix of coefficients
% b is nx1 column vector of constants
% x is nx1 column vector with the result
% Syntax: x = luab(A, b)

% Step 1: LU Factorization of A
[L,U] = lufactor(A);

[n,nc] = size(A);

% Step 2: Set Ly = b. Solve y using forward substitution
disp('-----------------------------')
disp('Forward Substitution')
disp('-----------------------------')
disp(' ')
% Initialize y
y = zeros(n,1); % initialize y as a column vector

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

% Step 3: Set Ux = y. Solve x using backward substitution
disp('-----------------------------')
disp('Backward Substitution')
disp('-----------------------------')
disp(' ')
% Initialize x
x = zeros(n,1); % initialize x as a column vector

x(n)=y(n)/U(n,n);
for i = n:-1:1
    sum = 0;
    for j = i+1:n
        sum = sum + U(i,j)*x(j);
    end
    x(i) = (y(i)-sum)/U(i,i);
end
disp(x)