> rk2:=proc(f1, f2, a, ya1, ya2,h,N)
> local t ,w1, w2, i, k11, k12, k13, k14, k21, k22, k23, k24;
> 
> t := a;
> w1[0] := ya1;
> w2[0] := ya2;
> 
> for i from 0 to N-1 do
>   k11 := evalf(h*f1(t, w1[i], w2[i]));
>   k21 := evalf(h*f2(t, w1[i], w2[i]));
> 
>   k12 := h*f1(t+h/2, w1[i]+k11/2, w2[i]+k21/2);
>   k22 := h*f2(t+h/2, w1[i]+k11/2, w2[i]+k21/2);
> 
>   k13 := h*f1(t+h/2, w1[i]+k12/2, w2[i]+k22/2);
>   k23 := h*f2(t+h/2, w1[i]+k12/2, w2[i]+k22/2);
> 
>   k14 := h*f1(t+h, w1[i]+k13, w2[i]+k23);
>   k24 := h*f2(t+h, w1[i]+k13, w2[i]+k23);
> 
>   w1[i+1] := w1[i] + (k11+2*(k12+k13)+k14)/6;
>   w2[i+1] := w2[i] + (k21+2*(k22+k23)+k24)/6;
> 
>   t := t + h
> end do;
> return(w1, w2)
> end proc;