#
#  Program that, using Runge-Kutta method, 
#  solves mxm system of IVP of ODE
#        x[i]' = f[i](t,x[1],x[2],...,x[m]),     t in [a,b]
#        x[i](0) = x0[i], 
#                i=1,2,...,m
#
#        input: 
#               f        --- a Vector of operators
#                            for example:
#                              > f:=Vector(2);
#                              > f[1]:=(t,x,y)->x-0.1*x*y; 
#                                f[2]:=(t,x,y)->-.5*y+.02*x*y;
#               a, b     --- end points of the t interval
#               x0       --- Vector of initial values
#               n        --- number of points  dividing [a,b]
#
#       output: t        --- the Vector of nodes dividing [a,b]
#               x        --- the array(1..m), 
#                            each entry x[i] is a Vector 
#                            x[i][j] is the approximate value 
#                            of x[i](t[j])
#
RungeKutta := proc( f :: Vector,   # vector of RHS functions
                    a :: numeric,  # left end of the interval
                    b :: numeric,  # right end of the interval
                   x0 :: Vector,   # vector of initial values
                    n :: integer   # number of points for t
                            )
   local i, j, k, l, m, delta, delta2, t, x;

   m := LinearAlgebra[Dimension](f);

   delta:=evalf( (b-a)/(n-1) );  
   delta2:=0.5*delta;

   t:=Vector([ seq( evalf( a+(i-1)*delta ), i=1..n )]);
   x:=array([ seq( Vector(n), i=1..m) ] );

   k:=array(1..4, [ seq( array(1..m), i=1..4) ] );

   for i to m do
      x[i][1]:=evalf(x0[i])
   end do;

   for i from 1 to n-1 do
      for j to m do
         k[1][j]:= evalf( delta*f[j](t[i], seq( x[l][i], l=1..m ) ) )
      end do;
      for j to m do
         k[2][j]:= evalf( delta*f[j](t[i]+delta2, 
                     seq( x[l][i]+0.5*k[1][l], l=1..m ) ) )
      end do;
      for j to m do
         k[3][j]:= evalf( delta*f[j](t[i]+delta2, 
                     seq( x[l][i]+0.5*k[2][l], l=1..m ) ) )
      end do;
      for j to m do
         k[4][j]:= evalf( delta*f[j](t[i]+delta,  
                     seq( x[l][i]+k[3][l], l=1..m ) ) )
      end do;
      for j to m do
         x[j][i+1]:= evalf( x[j][i] + 
                    (k[1][j]+2.0*(k[2][j]+k[3][j])+k[4][j])/6.0 )
      end do;
   end do;

   print(` End of Runge-Kutta method for mxm IVP-ODE `);

   return t, x;

end;