Code snippets for ODE Driver/Stepper

ODE23Step[y_,f_,h_]:= Module[
{k1,k2,k3,k4,yNew,zNew},
k1=f[y];
k2=f[y+ 1/2 h k1];
k3=f[y+3/4 h k2];
yNew=y + h ( 2/9 k1 + 1/3 k2 + 4/9 k3);
k4=f[yNew];
zNew=y +h*(7/24 k1 + 1/4 k2 + 1/3 k3 + 1/8 k4);
{yNew, Norm[zNew-yNew]}
]

ODE23StepBAD[y_,f_,h_]:= Module[
{k1,k2,k3,k4,yNew,zNew},
k1=f[y];
k2=f[y+ 1/2 h k1];
k3=f[y+3/4 h k2];
yNew=y + h ( 2/9 k1 + 1/3 k2 + 4/9 k3);
k4=f[yNew];
zNew=y +h*(7/24 k1 + 1/4 k2 + 1/3 k3 + 1/8 k4);
(* Returning the 2nd order solution *)
{zNew, Norm[zNew-yNew]}
]

Testing the thing it solve y'[t] =f[y[t]]

f[y_]:= Sin[(y^2)] + Cos[y] + {{-2, 1},{1,-1}}.y
y0={1.0,2.0};
h=0.1;
ConvData=Table[
h= (2.0^(-n));
{h,ODE23Step[y0,f,h][[2]]},
{n, 4,12}];
ListLogLogPlot[ConvData,
AxesLabel->{"h","error"}]
Clear[h]
Fit[Log[ConvData],{1,h},h]

So this is the correct behavior.

Building a driver for ODE23Step

Driver[y0_,f_,n_, TMax_,Stepper_]:=Module[
{Data,yy=y0,err,h=TMax/n,tt=0.0},
  Data = Table[
   {yy, err} = Stepper[yy, f, h];
   tt = tt + h;
   {tt, yy},
   {n}];
(* fixes it so that it comes out in the form {ts, ys} *)
Transpose[Data]
]

Clear[f,y,t,tt,yy,h,n,ySol]
f[y_]:= Sin[(y^2)] + Cos[y] + {{-0.5, 1},{-1,0}}.y
y0={1.0,2.0};
n=600;
TMax=220;
{tData,yData}=Driver[y0,f,n, TMax,ODE23Step];

ySol[t_]=y[t]/.NDSolve[ {
y'[t]==f[y[t]],
y[0]==y0
},
y[t],
{t, 0, TMax}][[1]];
MmaPic=ParametricPlot[ySol[t],{t, 0,TMax},
PlotRange->All,
AspectRatio->1];
OurPic=ListPlot[yData,
PlotRange->All,
AxesLabel->{"y1","y2"}
];
Show[OurPic, MmaPic]