(* RUNGE-KUTTA METHOD FOR SYSTEMS OF DIFFERENTIAL EQUATIONS ALGORITHM 5.7
*
*  To approximate the solution of the mth order system of 
*   first-order initial-value problems
*		UJ' = FJ( T, U1, U2, ..., UM), J = 1,2,3
*	  	A <= T <=B, UJ(A) = ALPHAJ, J = 1,2,3
*   at (n+1) equally spaced numbers in the interval [A,B].
*
*   INPUT:  endpoints a, b; number of equations M < 3; initial 
*            conditions ALPHA1, ..., ALPHA3; integer n.
*
*  OUTPUT:  approximation WJ to UJ(T) at the (n+1) values of T.
*)
OK = 0;
While[OK == 0,
   M = Input["Input the number of equations"];
   If[M <= 0 || M > 3,
      Input["Number must be a positive integer less than\n
      or equal to three/n"],
      OK = 1;
   ];
];
For[i = 1, i <= M, i++,
   TEMP[i]=Input["Input the function F("<>ToString[i]<>")(t,y1,y2,y3) in terms\n
   of t and y1,y2,y3\n
   \n
   for example: y1-t^2+1\n"];
   F[i][t_,y1_,y2_,y3_] :=Evaluate[TEMP[i]];
];
OK = 0;
While[OK == 0,
   A=Input["Input the left endpoint\n"];  
   B=Input["Input the right endpoint\n"];    
   If[A >= B,
      Input["Left endpoint must be less than right endpoint\n
      \n
      Press 1 [enter] to continue\n"],
      OK = 1;
   ];
];
For[i=0,
   i <= M-1,
   i=i+1;
   ALPHA[i]=Input["Input the initial condition alpha("<>ToString[i]<>")\n"];
];
OK=0;
While[OK == 0,
   n=Input["Input a positive integer for the number of\n
   subintervals\n"];
   If[n <= 0,
      Input["Number must be a positive integer\n
      \n
      Press 1 [enter] to continue\n"],
      OK = 1;
   ];
];
If[OK == 1,
   FLAG = Input["Select output destination\n
                 1. Screen\n
                 2. Text file\n
                 Enter 1 or 2\n"];
   If[FLAG == 2,
      NAME = InputString["Input the file name\n
                          For example:   output.dta\n"];
      OUP = OpenWrite[NAME,FormatType->OutputForm],
      OUP = "stdout";
   ];
   Write[OUP,"RUNGE-KUTTA METHOD FOR SYSTEMS\n
      OF DIFFERENTIAL EQUATIONS\n"];
   Write[OUP,"\n"];
   Write[OUP,"t   \n"];
   For[i=0,
      i <= M-1,
      i=i+1;
      Write[OUP,"W",i,"\n"]; 
   ];
   (* Step 1 *)
   H = (B-A)/n;
   T = A;
   (* Step 2 *)
   For[J=0,
      J <= M-1,
      J=J+1;
      W[J]=ALPHA[J];
   ];
   (* Step 3 *)
   Write[OUP,"\n"];
   Write[OUP,SetPrecision[T,9],"\n"];
   For[i=0,
      i <= M-1,
      i= i+1;
      Write[OUP,SetPrecision[W[i],9],"\n"];
   ];
   Write[OUP,"\n"];
   (* Step 4 *)
   For[L=0,
      L <= n-1,
      L=L+1;
      (* Step 5 *)
      For[J=0,
	 J <= M-1,
	 J = J+1;
         K[1,J]=N[H*F[J][T,W[1],W[2],W[3]]];
      ];
      (* Step 6 *)
      For[J=0,
	 J <= M-1,
	 J = J+1;
         K[2,J]=N[H*F[J][T+H/2.0,W[1]+K[1,1]/2.0,W[2]+K[1,2]/2.0,W[3]+K[1,3]/2.0]];
      ];
      (* Step 7 *)
      For[J=0,
	 J <= M-1,
	 J = J+1;
         K[3,J]=N[H*F[J][T+H/2.0,W[1]+K[2,1]/2.0,W[2]+K[2,2]/2.0,W[3]+K[1,3]/2.0]];
      ];
      (* Step 8 *)
      For[J=0,
	 J <= M-1,
	 J = J+1;
         K[4,J]=N[H*F[J][T+H,W[1]+K[3,1],W[2]+K[3,2],W[3]+K[3,3]]];
      ];
      (* Step 9 *)
      For[J=0,
	 J <= M-1,
	 J = J+1;
         W[J]=W[J]+(K[1,J]+2.0*K[2,J]+2.0*K[3,J]+K[4,J])/6.0;
      ];
      (* Step 10 *)
      T=A+L*H;
      (* Step 11 *)
      Write[OUP,SetPrecision[T,9],"\n"];
      For[i=0,
         i <= M-1,
         i= i+1;
         Write[OUP,SetPrecision[W[i],9],"\n"];
      ];
      Write[OUP,"\n"];
   ];
   (* Step 12 *)
   If[OUP == "OutputStream[",NAME," 3]",
      Print["Output file: ",NAME," created successfully\n"];
      Close[OUP];
   ];
];
Quit[];