>  restart; 
> with(LinearAlgebra):
> #Vorwaertselimination, ueberschreibt b
>  velim:=proc(Umat,b) 
>  local j,k,n; 
>  n:=Dimension(b); 
>  for j from 1 to n do 
>     for k from 1 by 1 to j-1 do b[j]:=b[j] - Umat[j,k]*b[k] end do; 
>     b[j]:=b[j]/Umat[j,j];
>  end do; 
>  print(b);
>  end;
>  L:=Matrix([[1,0,0,0],[2,3,0,0],[4,5,6,0],[7,8,9,10]]);
>  b:=Vector([1,2,3,4]);
>  velim(L,b);
>                     
>