clc
clear all
close all

function x = gauss_elimination(A, b)

     % Solve A*x = b using Gaussian elimination with partial pivoting
     [n, m] = size(A);
     if n ~= m
          error("Matrix A must be square");
     end
     b = b(:);           % ensure column vector
     if length(b) ~= n
          error("Size of b must match A");
     end

     % Forward elimination
     for k = 1:n-1
          % Pivoting: find row with largest |A(i,k)|, i >= k
          [~, piv] = max(abs(A(k:n, k)));
          piv = piv + k - 1;
          if A(piv, k) == 0
               error("Matrix is singular or nearly singular");
          end
          % Swap rows in A and b if needed
          if piv ~= k
               A([k, piv], :) = A([piv, k], :);
               b([k, piv])   = b([piv, k]);
          end

          % Eliminate entries below pivot
          for i = k+1:n
               factor = A(i, k) / A(k, k);
               A(i, k:m) = A(i, k:m) - factor * A(k, k:m);
               b(i)      = b(i)      - factor * b(k);
          end
     end

     % Back substitution
     x = zeros(n, 1);
     for i = n:-1:1
          s = A(i, i+1:n) * x(i+1:n);
          x(i) = (b(i) - s) / A(i, i);
     end
end

gauss_elimination([1 2 3; 4 7 8; 1 1 2], [1 4 6]')

[1 2 3; 4 7 8; 1 1 2]\[1; 4; 6]