clear all;

n = 4;

a = [16 -8 -4 0;
    -8 13 -4 -3;
    -4 -4 9 0;
    0 -3 0 3];

u = [0 0 0 0;
    0 0 0 0;
    0 0 0 0;
    0 0 0 0];

for i = 1:1:n

  s = 0;
  for k = 1:1:(i-1)
    s += u(k,i)^2;
  endfor
  u(i,i) = sqrt(a(i,iu) - s);

  for j = (i+1):1:n

    s = 0;
    for k = 1:1:(i-1)
      s += u(k,i)*u(k,j);
    endfor

    u(i,j) = (a(i,j) - s) / u(i,i);
  endfor
endfor