function p3
  A = [1,4,3;2,7,9;5,8,-2];
  gausslu(A);
end
% -----------------------------------------------------------
function [L, U] = gausslu(A)
  U = A;
  
  [n, dummy] = size(A);
  L = eye(n);

  for step = 1:n-1
    pivot = U(step, step);
    for i = step+1:n
      factor = U(i,step)/pivot;
      L(i,step) = factor;
      U(i,step+1:n) = U(i,step+1:n) - factor*U(step,step+1:n);
      U(i,step) = 0;
    end
  end

  L
  U
  LU = L*U
  A
end
