function A=LU_golub_kij(A)
[m,n]=size(A);
A = sym(A);
disp('Begin');
for k=1:min(m, n)-1
    for i=k+1:m
        A(i,k)=A(i,k)/A(k,k);
        for j=k+1:n
           A(i, j) = A(i, j) - A(i,k)*A(k, j); 
        end
    end
end    
disp(A);
disp('End');