restart;
 
with(LinearAlgebra):
 
loesAxb:=proc(A::Matrix,b::Vector)

 local i,n,gaussTp,vorelim1,ruecksubs,x;
 n:=Dimension(b);
 x:=Vector(n);
 gaussTp:= proc(A)
 global ps,pz;
 local j,k,m,n,piv,pivj,pivk,x;
 n:=RowDimension(A);
 x:=Vector(n);
 ps:=Vector(n,i->i);
 pz:=Vector(n,i->i);
 for m from 1 to n do
     # Pivotsuche
     piv:=abs(A[m,m]);
     pivj:=m;pivk:=m;
     for k from m to n do
         for j from m to n do
             if abs(A[j,k])>piv then
                piv:=abs(A[j,k]);
                pivj:=j;pivk:=k
             end if
         end do
     end do;
     # Zeilentausch
     j:=pz[m]; pz[m]:=pz[pivj]; pz[pivj]:=j;
     for k from 1 to n do
         x[k]:=A[m,k]; A[m,k]:=A[pivj,k];
         A[pivj,k]:=x[k]
     end do;
     # Spaltentausch
     k:=ps[m]; ps[m]:=ps[pivk]; ps[pivk]:=k;
     for j from 1 to n do
         x[j]:=A[j,m]; A[j,m]:=A[j,pivk];
         A[j,pivk]:=x[j]
     end do;
     # Eliminationsschritt
     for j from m+1 to n do
        A[j,m]:=A[j,m]/A[m,m];
        for k from m+1 to n do
           A[j,k]:=A[j,k]-A[j,m]*A[m,k]
        end do
     end do
 end do
 end proc;
  
 vorelim1:=proc(L,b)
  local j,k,n;
  n:=Dimension(b);
  for j from 1 to n do
     for k from 1 to j-1 do
         b[j]:=b[j]-L[j,k]*b[k]
     end do
   end do 
 end proc;
  
 ruecksubs:=proc(U,b)
  local j,k,n;
  n:=Dimension(b);
  for j from n by (-1) to 1 do 
     for k from n by (-1) to j+1 do
         b[j]:=b[j]-U[j,k]*b[k]
     end do;
     b[j]:=b[j]/U[j,j] 
  end do
 end proc;
  
 gaussTp(A);
 # b-Zeilentausch
   for i from 1 to n do
      x[i]:=b[pz[i]]
   end do;
 vorelim1(A,x);
 ruecksubs(A,x);
 # Spaltentausch rueckgaengig
   for i from 1 to n do
       b[ps[i]]:=x[i]
   end do;
 print(b)
end proc:

> A:=Matrix([[1,1],[99999999,100000000]]);
  b:=Vector([1,9]);
  loesAxb(A,b);