RK4[f_,u0_,h_,t0_,T_]:=(
n=Ceiling[(T-t0)/h];
t=Table[i h,{i,0,n}];
y=Table[0,{n+1}];
y[[1]]=u0;
For[i=1,i≤n,i++,
k1=h f[t[[i]],y[[i]]];
k2=h f[t[[i]]+h/2,y[[i]]+k1/2];
k3=h f[t[[i]]+h/2,y[[i]]+k2/2];
k4=h f[t[[i]]+h,y[[i]]+k3];
y[[i+1]]=y[[i]]+(1/6)(k1+2k2+2k3+k4)
];
y
)