> restart; 
> with(LinearAlgebra): 
> printlevel:=0:
> 
> #Jacobi-Verfahren , Maxnorm des Residuums soll kleiner tol sein
> Jacobi:=proc(A,b,tol) 
> global res,vnorm,x ;
> local n,k, vectornorm, residuum, JacobiIt;
> # option trace;
> JacobiIt:=proc(A,x,b) #Jacobi-Iterationsschritt
> local j,k,n,y,z;
> # option trace;
> n:=Dimension(b); 
> y:=Vector(n);
> for j from 1 to n do
>   z:=b[j];
>   for k from 1 to j-1 do z:=z - A[j,k]*x[k] end do;
>   for k from j+1 to n do z:=z - A[j,k]*x[k] end do;
>   y[j]:=z/A[j,j]
> end do;
> for j from 1 to n do x[j]:=y[j] end do;
> end proc: 
> vectornorm:=proc(x)  #Maximumnorm des Vektors x
> global vnorm;
> local j, n;
> # option trace;
> vnorm:=0; n:=Dimension(x); 
> for j from 1 to n do
>  if abs(x[j]) >vnorm then vnorm :=abs(x[j]) end if
> end do;
> #eval(vnorm)
> end proc;
> residuum:=proc(A,x,b)   #Residuum Ax-b
> global res;
> local j,k,n;
> # option trace;
> n:=Dimension(b); res:=Vector(n);
> for j from 1 to n do
>   for k from 1 to n do res[j]:= res[j]+A[j,k]*x[k] end do;
>   res[j]:=res[j]-b[j]
> end do;
> #evalm(res)
> end proc;
> n:=Dimension(b); x:=Vector(n); res:=Vector(n); 
> residuum(A,x,b);
> vnorm:= vectornorm(res); k:=0;
> while (vnorm > tol) do
>   JacobiIt(A,x,b);
>   residuum(A,x,b);
>   vectornorm(res); 
>   k:=k+1;
>   print(k, x); # druckt nach jedem Iterationsschritt Naeherung x
> end do; 
> end proc:
> A:=Matrix([[10.,-1.,0.,2.],[1.,12.,-1.,2.],[-2.,1.,15.,0.],[1.,-2.,0.,20.]]):
> b:=Vector([11,14,14,19]): tol:=0.001:  Jacobi(A,b,tol);
>                                        
>