function x = lusolve(L, U, b)
%LUSOLVE - solve the system of linear equations A*x = b
%          using LU-decomposition

% This version of the code uses vectorization

    [n, ~] = size(L);
    y = zeros(n, 1);

    % Solve L y = b
    for k = 1:n
        x(k) = b(k);
        i = k+1:n;
        b(i) = b(i) - L(i,k)*x(k);
    end

    % preallocate the solution vector
    x = zeros(n, 1);

    % Solve u*x = y
    for k = n:-1:1
        x(k) = b(k)/U(k,k);
        i = 1:k-1;
        b(i) = b(i) - U(i,k)*x(k);
    end

end
