function [L,U] = mylu1(A)
% LU factorization

n = length(A);
for k=1:n-1
    if A(k,k) == 0, error('zero division'), end
    % part of L
    for i=k+1:n
        A(i,k) = A(i,k)/A(k,k);
    end
    % part of U
    for j=k+1:n
        for i=k+1:n
            A(i,j) = A(i,j) - A(i,k)*A(k,j);
        end
    end
end
L = eye(n) + tril(A,-1);
U = triu(A);