% declaration of functions

function a = gepivot(a, k)
[n, m] = size(a);
piv = a(k,k);
newpivot = k;

for i = k+1:n            % check all rows below pivot
   if abs(a(i,k)) > abs(piv)
       newpivot = i;
       piv = a(i,k);
   end
end

if piv == 0              % trouble - can't avoid 0 pivot
   error('Matrix is singular');
else
 if newpivot ~= k     % best pivot elsewhere, swap rows
   for j = k:m
     temp1 = a(k,j);
     temp2 = a(newpivot,j);
     a(k,j) = temp2;
     a(newpivot,j) = temp1;
   end
end
end
end

function x = gauselim(a, b)
% Gaussian Elimination routine
%
% USAGE:   x = gauselim(a, b)
%  where   a = square coefficient matrix, size n X n
%          b = single RHS column vector, size n X 1
%          x = solution column vector, size n X 1


% check that input is legal
[n,m] = size(a);
if n ~= m
  error('Error in gauselim: Matrix must be square.');
 end

%set up space for solution vector
[nlign,ncol]=size(b);
x = zeros(nlign,ncol);

% form augmented matrix by adding b to last column of a
m = n+1;
a(:,m) = b;

% FORWARD ELIMINATION
for k = 1:n-1

  a = gepivot(a, k);   % partial pivoting

  
  for i = k+1:n
    term = a(i,k)/a(k,k);
    for j = k+1:m
       a(i,j) = a(i,j) - term*a(k,j);
    end
  end
end

% BACK SUBSTITUTION
if a(n,n) == 0
     error('Matrix is singular');
   end

x(n) = a(n,m)/a(n,n);
for i = n-1:-1:1
   sum = 0;
   for j = i+1:n
      sum = sum + a(i,j)*x(j);
   end
   x(i) = (a(i,m) - sum) / a(i,i);
end
end

% start of the program

% input matrix a
a = [1,  2;
     3,  4];

% input vector b
b = [5, 6];

% transpose of vector b
t=b';

% call to the gauss elimination subroutine
ans = gauselim(a, t);

% display the result
write("answer = ", ans);