[MATLAB] Solving Ax = b with Factorization A = LU

function [L, U] = factorize(A)
% square LU factorization with no row exchanges
[n, n] = size(A); 
zr = 1.e-6;

for k = 1:n
    if abs(A(k,k)) < zr
    end 
    L(k, k) = 1;   %set 1s on diagonal 
    for i = k+1:n
        L(i, k) = A(i, k)/A(k, k);   % computes multipliers
        for j = k+1:n
            A(i, j) = A(i, j) - L(i, k)*A(k, j);   % substracts mulitipliers times A(k) from A(i)
        end
    end
    for j = k:n
        U(k, j) = A(k, j);    % upper triangular matrix 
    end
end 

function x = slv(A, b)
% Solve Ax = b using Lc=b and Ux=c
[L, U] = factorize(A)
[n, n] = size(A);
s = 0;  
t = 0;

for k=1:n     % Forward Elimination Lc = b to slove c 
    for j=1:k-1
        s = s + L(k, j)*c(j);
    end
    c(k) = b(k) - s;
    s=0;
end

for k=n:-1:1     % Back-Substitution with Ux = c to solve from x(n) to x(1)
    for j=k+1:n     % from diagnoal to end of row
        t = t + U(k, j)*x(j);  % adds up U times previously earned x 
    end
    x(k) = (c(k) - t)/U(k, k);  % divide by pivot 
end
x = x';