> restart;
> with(LinearAlgebra):
> 

> Inv:=proc(Umat::Matrix) 
> local j,r,k,n,M, summe; 
> n:=RowDimension(Umat); 
> r:=0; 
> for j from 1 by 1 to n do 
>    if Umat[j,j]=0 then  r:=1 end if 
> end do; 
> if r=1 then print(`Matrix nicht invertierbar`) 
> else M:=matrix(n,n); 
>      for j from 1 by 1 to n do 
>        for k from 1 by 1 to n do M[j,k]:=0 end do 
>      end do; 
>      for r from 1 by 1 to n do  
>         M[r,r]:=1/Umat[r,r]; 
>         for j from r+1 by 1 to n do 
>            summe:=0; 
>            for k from r by 1 to j-1 do summe:= summe + Umat[j,k]*M[k,r] end do; 
>            M[j,r]:=-summe/Umat[j,j] 
>         end do
>      end do; 
>      print(`Inverse=`, M) 
> end if 
> end proc:
> L:=Matrix([[1,0,0,0],[2,3,0,0],[4,5,6,0],[7,8,9,10]]); Inv(L);