> #Einschrittverfahren
> restart;
> with(linalg):
> 
> #vektorielle Version
> RungeKuttasystem:=proc(y0,t0,f,T,N)
> # Schrittweite (T-t0)/N
> # option trace;
> local d,k,y,t,k1,k2,k3,k4,j,h;
> d:=vectdim(y0);
> y:=vector(d); 
> k1:=vector(d); k2:=vector(d); k3:=vector(d); k4:=vector(d);
> for j from 1 to d do y[j]:=y0[j] od; 
> t:=t0; h:=(T-t0)/N;
> for j from 0 to N-1 do
>  k1:=f(t,y); 
>  k2:=f(t+h/2,matadd(y,k1,1,h/2)); 
>  k3:=f(t+h/2,matadd(y,k2,1,h/2));
>  t:=t+h; 
>  k4:=f(t,matadd(y,k3,1,h)); 
>  for k from 1 to d do y[k]:=y[k]+h*(k1[k]+2*k2[k]+2*k3[k]+k4[k])/6 od;
> print(t,y)
> od
> end:
> 
Warning, the protected names norm and trace have been redefined and unprotected

> #Beispiel:
> f1:=proc(t,y) local g;
> g:=vector(2);
> g[1]:=y[2]; g[2]:=2*y[1]^3;
> evalm(g) end:
> 
> y0:=vector(2,[1.0,-1.0]):
> t0:=-2.0:
> T:=-1.0:
> N:=10:
> RungeKuttasystem(y0,t0,f1,T,N);
> # exakte Lösung:
> for j from 1 to 10 do print(-2+j*0.1,evalf(1/(-2+j*0.1+3))) od;
> 

              -1.900000000, [.9090945296, -.8264416310]


              -1.800000000, [.8333395758, -.6944356321]


              -1.700000000, [.7692392619, -.5917032693]


              -1.600000000, [.7142964020, -.5101876015]


              -1.500000000, [.6666796664, -.4444242155]


              -1.400000000, [.6250155298, -.3906009810]


              -1.300000000, [.5882536367, -.3459928660]


              -1.200000000, [.5555770389, -.3086100862]


              -1.100000000, [.5263407762, -.2769722875]


              -1.000000000, [.5000288815, -.2499596877]


                          -1.9, .9090909091


                          -1.8, .8333333333


                          -1.7, .7692307692


                          -1.6, .7142857143


                          -1.5, .6666666667


                          -1.4, .6250000000


                          -1.3, .5882352941


                          -1.2, .5555555556


                          -1.1, .5263157895


                          -1.0, .5000000000

> 
>