> 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;
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 + 1/2*h, w1[i] + 1/2*k11, w2[i] + 1/2*k21);
    k22 := h*f2(t + 1/2*h, w1[i] + 1/2*k11, w2[i] + 1/2*k21);
    k13 := h*f1(t + 1/2*h, w1[i] + 1/2*k12, w2[i] + 1/2*k22);
    k23 := h*f2(t + 1/2*h, w1[i] + 1/2*k12, w2[i] + 1/2*k22);
    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] + 1/6*k11 + 1/3*k12 + 1/3*k13 + 1/6*k14;
    w2[i + 1] := w2[i] + 1/6*k21 + 1/3*k22 + 1/3*k23 + 1/6*k24;
    t := t + h;
  end do;
  return w1, w2;
end proc;