> #runge-kutta for ode
> myrk:=proc(f,t0,y0,tstep,n)
> local told, yold, tnew, ynew, k1, k2, k3, k4, deltay, yprime, i;
> tnew:=t0;ynew:=y0;
> printf("t, y, k1, k2, k3, k4, y', delta y, new y");
> for i from 1 to n do
> told:=tnew;yold:=ynew;
> k1:=eval(f,{t=told,y=yold});
> k2:=eval(f,{t=told+tstep/2,y=yold+tstep*k1/2});
> k3:=eval(f,{t=told+tstep/2,y=yold+tstep*k2/2});
> k4:=eval(f,{t=tnew,y=yold+tstep*k3});
> yprime:=(k1+2*k2+2*k3+k4)/6;
> deltay:= tstep*yprime;
> tnew:=told+tstep;ynew:=yold+deltay;
> print(told,yold,k1,k2,k3,k4,yprime,deltay,ynew);
> od
> end:
> myrk(1+y^2,0,0,0.5,5);
t, y, k1, k2, k3, k4, y', delta y, new y

  0, 0, 1, 1.062500000, 1.070556641, 1.286522880, 1.092106027,

        0.5460530135, 0.5460530135


  0.5, 0.5460530135, 1.298173894, 1.757938243, 1.971284310,

        3.346090088, 2.017118181, 1.008559090, 1.554612104


  1.0, 1.554612104, 3.416818794, 6.802398386, 11.59640321,

        55.06386944, 15.87971524, 7.939857620, 9.494469724


  1.5, 9.494469724, 91.14495534, 1043.044145, 73039.03847,

                       10                 9                 9
        0.1334368844 10  , 0.2224195165 10 , 0.1112097582 10 ,

                       9
        0.1112097677 10


                      9                 17                 31
  2.0, 0.1112097677 10 , 0.1236761243 10  , 0.9559865516 10  ,

                       61                 121                 121
        0.5711939293 10  , 0.8156562619 10   , 0.1359427103 10   ,

                       120                 120
        0.6797135515 10   , 0.6797135515 10

> myrk([y[2],-y[1]],0,[1,-1],0.1,3);# This is y''-y=0, y(0)=1, y'(0)=-1 translated to system
t, y, k1, k2, k3, k4, y', delta y, new y

  0, [1, -1], [-1, -1], [-1.050000000, -0.9500000000],

        [-1.047500000, -0.9475000000], [-1.094750000, -0.8952500000],

        [-1.048291667, -0.9483750000],

        [-0.1048291667, -0.09483750000], [0.8951708333, -1.094837500]


  0.1, [0.8951708333, -1.094837500], [-1.094837500, -0.8951708333],

        [-1.139596042, -0.8404289583], [-1.136858948, -0.8381910312],

        [-1.178656603, -0.7814849385], [-1.137734014, -0.8389826252],

        [-0.1137734014, -0.08389826252], [0.7813974319, -1.178735763]


  0.2, [0.7813974319, -1.178735763], [-1.178735763, -0.7813974319],

        [-1.217805635, -0.7224606438], [-1.214858795, -0.7205071502],

        [-1.250786478, -0.6599115524], [-1.215808517, -0.7212074287],

        [-0.1215808517, -0.07212074287], [0.6598165802, -1.250856506]

>