> restart; with(LinearAlgebra):
> # Loesen eines Gleichungssystems (ohne Pivot);
> loesAxb:=proc(A,b)
> local gauss1, vorelim1, ruecksubs;
> gauss1:=proc(A::Matrix) 
> #Gauss-Elimination mit Ueberschreiben von A, ohne Pivot;
> local j,k,m,n;
> # option trace;
> n:=RowDimension(A); 
> for m from 1 by 1 to n do 
>    for j from m+1 by 1 to n do 
>       A[j,m]:=A[j,m]/A[m,m]; 
>       for k from m+1 by 1 to n do A[j,k]:=A[j,k]- A[j,m]*A[m,k] end do; 
>    end do 
> end do;
> evalm(A)
> end proc; 
> # Vorwaertselimination, ueberschreibt b
> # Annahme: Diagonale 1 !!;
> vorelim1:=proc(L::Matrix,b::Vector) 
> local n,j,k; 
> # option trace;
> n:=Dimension(b); 
> for j from 1 by 1 to n do 
>    for k from 1 by 1 to j-1 do b[j]:=b[j]-L[j,k]*b[k] end do 
> end do; 
> evalm(b)
> end proc;
> ruecksubs:=proc(U::Matrix,b::Vector) 
> local n, j, k; 
> # option trace
> 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;
> evalm(b)
> end proc;
> gauss1(A); vorelim1(A,b);
> ruecksubs(A,b); 
> evalm(b);
> end proc:
> 
> A:=Matrix([[6.,7.,0,1,2,5],[13,4,2,4,5,3],[6,2,3,6,8,9],[1,3,0,0,2,4],[1,1,3,5,7,9],[0,0,4,5,0,2]]);
> b:=Vector([1,0,1,0,1,0]);
> loesAxb(A,b);
> #Kontrolle
> A:=Matrix([[6.,7.,0,1,2,5],[13,4,2,4,5,3],[6,2,3,6,8,9],[1,3,0,0,2,4],[1,1,3,5,7,9],[0,0,4,5,0,2]]);
> b:=Vector([1,0,1,0,1,0]);
> LinearSolve(A,b);
>