> restart;
> # RUNGE-KUTTA 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, ..., M
> #            A <= T <= B, UJ(A) = ALPHAJ, J = 1, 2, ..., M
> # AT (N+1) EQUALLY SPACED NUMBERS IN THE INTERVAL [A,B].
> #
> # INPUT:   ENDPOINTS A,B; NUMBER OF EQUATIONS M; INITIAL
> #          CONDITIONS ALPHA1, ..., ALPHAM; INTEGER N.
> #
> # OUTPUT:  APPROXIMATION WJ TO UJ(T) AT THE (N+1) VALUES OF T.
> alg057 := proc() local OK, M, I, F, A, B, ALPHA, N, FLAG, NAME, OUP, H, T, J, W, L, K;
> printf(`This is the Runge-Kutta Method for Systems of m equations\n`);
> OK := FALSE;
> while OK = FALSE do
> printf(`Input the number of equations\n`);
> M := scanf(`%d`)[1];
> if M <= 0 then 
> printf(`Number must be a positive integer\n`);
> else
> OK := TRUE;
> fi;
> od;
> for I from 1 to M do
> printf(`Input the function F[%d](t,y1 ... y%d) in terms of t\n`,I,M); 
> printf(`and y1 ... y%d\n`, M);
> printf(`For example: y1-t^2+1\n`);
> F[I] := scanf(`%a`)[1];
> F[I] := unapply(F[I],t,evaln(y || (1..M)));
> od;
> OK := FALSE;
> while OK = FALSE do
> printf(`Input left and right endpoints separated by blank\n`);
> A := scanf(`%f`)[1];
> B := scanf(`%f`)[1];
> if A >= B then 
> printf(`Left endpoint must be less than right endpoint\n`);
> else
> OK := TRUE;
> fi;
> od;
> for I from 1 to M do
> printf(`Input the initial condition alpha[%d]\n`, I);
> ALPHA[I] := scanf(`%f`)[1];
> od;
> OK := FALSE;
> while OK = FALSE do
> printf(`Input a positive integer for the number of subintervals\n`);
> N := scanf(`%d`)[1];
> if N <= 0 then
> printf(`Number must be a positive integer\n`);
> else
> OK := TRUE;
> fi;
> od;
> if OK = TRUE then
> printf(`Choice of output method:\n`);
> printf(`1. Output to screen\n`);
> printf(`2. Output to text file\n`);
> printf(`Please enter 1 or 2\n`);
> FLAG := scanf(`%d`)[1];
> if FLAG = 2 then
> printf(`Input the file name in the form - drive:\\name.ext\n`);
> printf(`For example   A:\\OUTPUT.DTA\n`);
> NAME := scanf(`%s`)[1];
> OUP := fopen(NAME,WRITE,TEXT);
> else
> OUP := default;
> fi;
> fprintf(OUP, `RUNGE-KUTTA METHOD FOR SYSTEMS OF DIFFERENTIAL EQUATIONS\n\n`);
> fprintf(OUP, `    T`);
> for I from 1 to M do
> fprintf(OUP, `          W%d`, I);
> od;
> # Step 1
> H := (B-A)/N;
> T := A;
> # Step 2
> for J from 1 to M do
> W[J] := ALPHA[J];
> od;
> # Step 3
> fprintf(OUP, `\n%5.3f`, T);
> for I from 1 to M do
> fprintf(OUP, ` %11.8f`, W[I]);
> od;
> fprintf(OUP, `\n`);
> # Step 4
> for L from 1 to N do
> # Step 5
> for J from 1 to M do
> K[1,J] := H*F[J](T,seq(W[i],i=1..M));
> od;
> # Step 6
> for J from 1 to M do
> K[2,J] := H*F[J](T+H/2.0,seq(W[i]+K[1,i]/2.0,i=1..M));
> od;
> # Step 7
> for J from 1 to M do
> K[3,J] := H*F[J](T+H/2.0,seq(W[i]+K[2,i]/2.0,i=1..M));
> od;
> # Step 8
> for J from 1 to M do
> K[4,J] := H*F[J](T+H,seq(W[i]+K[3,i],i=1..M));
> od;
> # Step 9
> for J from 1 to M do
> W[J] := W[J]+(K[1,J]+2.0*K[2,J]+2.0*K[3,J]+K[4,J])/6.0;
> od;
> # Step 10
> T := A+L*H;
> # Step 11
> fprintf(OUP, `%5.3f`, T);
> for I from 1 to M do
> fprintf(OUP, ` %11.8f`, W[I]);
> od;
> fprintf(OUP, `\n`);
> od;
> # Step 12
> if OUP <> default then
> fclose(OUP):
> printf(`Output file %s created successfully`,NAME);
> fi;
> fi;
> RETURN(0);
> end;
> alg057();