% gauss.m: Gaussian elimination (without pivoting) and backward substitution. 

  clear;
  format short;

% Step 1: Define the matrix A and the vector b
 
  n = 4;
  A = [4, 3, 2, 1;
       3, 4, 3, 2;
       2, 3, 4, 3;
       1, 2, 3, 4];
  b = [1, 1, -1, -1]; 
  
% Step 2: Factorization. 

  L = zeros(n,n);   % initialize the matrix L
  U = zeros(n,n);   % initialize the matrix U
  for k = 1:n
    L(k,k) = 1;     % compute L(k,k)
  end  
  for k = 1:n-1
    for i = k+1:n
      (fill out one line to compute multipliers L(i,k). 
       Don't forget the semicolon.)
      for j = k+1:n
        (fill out one line to update A(i,j). 
         Don't forget the semicolon.)
      end 
      (fill out one line to update b(i). 
       Don't forget the semicolon.)
    end
  end 

  for i = 1:n
    for j = i:n
      U(i,j) = A(i,j);    % compute U(i,j)
    end 
  end 

% Step 3: back substitute. 

  x = zeros(n,1);         % initialize x
  x(n) = b(n)/U(n,n);     % compute x(n)
  for k = n-1:-1:1
    tmp = 0;              % tmp is a temporal variable
    for j = k+1:n
      tmp = tmp + U(k,j)*x(j);
    end
    x(k) = (b(k)-tmp)/U(k,k);  % compute other x(k)
  end 

% Step 4: Output solutions. 

  L                 % output L
  U                 % output U
  x                 % output x