> restart:
> # CONTINUATION METHOD FOR SYSTEMS ALGORITHM 104
> #
> # To approximate the solution of the nonlinear system F(X)=0 given
> # an initial approximation X:
> #
> # INPUT:   Number n of equations and unknowns: initial approximation
> #          X=(X(1),...,X(n)): number of Runge-Kutta 4 iterations N.
> #
> # OUTPUT:  Approximate solution X=(X(1),...,X(n)). 
> #          
> LINSYS := proc(N,OK,A,Y) local K, I2, Z, IR, IA, J, C, L, JA:
> K := N-1:
> OK := TRUE:
> I2 := 1:
> while OK = TRUE and I2 <= K do
> Z := abs(A[I2-1,I2-1]):
> IR := I2:
> IA := I2+1:
> for J from IA to N do
> if abs(A[J-1,I2-1]) > Z then
> IR := J:
> Z := abs(A[J-1,I2-1]):
> fi:
> od:
> if Z <= 1.0e-20 then
> OK := FALSE:
> else
> if IR <> I2 then
> for J from I2 to N+1 do
> C := A[I2-1,J-1]:
> A[I2-1,J-1] := A[IR-1,J-1]:
> A[IR-1,J-1] := C:
> od:
> fi:
> for J from IA to N do
> C :=A[J-1,I2-1]/A[I2-1,I2-1]:
> if abs(C) <= 1.0e-20 then
> C := 0:
> fi:
> for L from I2 to N+1 do
> A[J-1,L-1] := A[J-1,L-1]-C*A[I2-1,L-1]:
> od:
> od:
> fi:
> I2 := I2+1:
> od:
> if OK = TRUE
> then if abs(A[N-1,N-1]) <= 1.0e-20 then
> OK := FALSE:
> else
> Y[N-1] := A[N-1,N]/A[N-1,N-1]:
> for I2 from 1 to K do
> J := N-I2:
> JA := J+1:
> C := A[J-1,N]:
> for L from JA to N do
> C := C-A[J-1,L-1]*Y[L-1]:
> od:
> Y[J-1] := C/A[J-1,J-1]:
> od:
> fi:
> fi:
> if OK = FALSE then
> print(`Linear system is singular`):
> fi:
> end:
> print(`This is the Continuation Method for Nonlinear Systems.`):
> OK := FALSE:
> while OK = FALSE do
> print(`Input the number n of equations.`):
> N := scanf(`%d`)[1]: print(`n = `):print(N):
> if N >= 2 then
> OK := TRUE:
> else
> print(`n must be an integer greater than 1.\n`):
> fi:
> od:
> for I2 from 1 to N do
> print(`Input the function F_I in terms of x1 ... xn for I = `,I2):
> F[I2] := scanf(`%a`)[1]:  print(`Function is `): print(F[I2]):
> od:
> for I2 from 1 to N do
> for J from 1 to N do
> P[I2,J] := diff(F[I2],evaln(x || J)):
> P[I2,J] := unapply(P[I2,J],evaln(x || (1..N))):
> od:
> od:
> for I2 from 1 to N do
> F[I2] := unapply(F[I2],evaln(x || (1..N))):
> od:
> OK := FALSE:
> while OK = FALSE do
> print(`Input the number N for RK4.`):
> NN := scanf(`%d`)[1]: print(`N = `):print(NN):
> if NN > 0 then
> OK := TRUE:
> else
> print(`Must be a positive integer.`):
> fi:
> od:
> for I2 from 1 to N do
> print(`Input initial approximation X(I) for I = `, I2):
> X[I2-1] := scanf(`%f`)[1]: print(`Input is `): print(X[I2-1]):
> od:
> if OK = TRUE then
> print(`Select output destination`):
> print(`1. Screen`):
> print(`2. Text file`):
> print(`Enter 1 or 2`):
> FLAG := scanf(`%d`)[1]: print(`Input is `):print(FLAG):
> if FLAG = 2 then
> print(`Input the file name in the form - drive\\:name.ext`):
> print(`for example   A:\\OUTPUT.DTA`):
> NAME := scanf(`%s`)[1]: print(`Output file is `):print(NAME):
> OUP := fopen(NAME,WRITE,TEXT):
> else
> OUP := default:
> fi:
> fprintf(OUP,`This is the Continuation Method using the Runge Kutta Method of Order 4 \n\n`):
> fprintf(OUP,`N = %2d \n\n`,NN):
> fprintf(OUP,`Iter#  Approximation \n\n`):
> K:=0:
> fprintf(OUP, ` %2d`, K):
> for I2 from 1 to N do
> fprintf(OUP, ` %11.8f `, X[I2-1]):
> od:
> fprintf(OUP,` \n`):
> # Step 1
> H := 1/NN:
> for I2 from 1 to N do
> b[I2-1] := H*evalf(-F[I2](seq(X[i],i=0..N-1))):
> od:
> # Step 2
> for K from 1 to NN do
> # Steps 3 - 6
> for I2 from 1 to N do
> for J from 1 to N do
> A[I2-1,J-1] := evalf(P[I2,J](seq(X[i],i=0..N-1))):
> od:
> od:
> for I2 from 1 to N do
> A[I2-1,N] := b[I2-1]:
> od:
> LINSYS(N,OK,A,Y):
> if OK = FALSE then
> break:
> fi:
> for I2 from 1 to N do
> K1[I2-1] := Y[I2-1]:
> X1[I2-1] := X[I2-1]+0.5*K1[I2-1]:
> od:
> for I2 from 1 to N do
> for J from 1 to N do
> A[I2-1,J-1] := evalf(P[I2,J](seq(X1[i],i=0..N-1))):
> od:
> od:
> for I2 from 1 to N do
> A[I2-1,N] := b[I2-1]:
> od:
> LINSYS(N,OK,A,Y):
> if OK = FALSE then
> break:
> fi:
> for I2 from 1 to N do
> K2[I2-1] := Y[I2-1]:
> X2[I2-1] := X[I2-1]+0.5*K2[I2-1]:
> od:
> for I2 from 1 to N do
> for J from 1 to N do
> A[I2-1,J-1] := evalf(P[I2,J](seq(X2[i],i=0..N-1))):
> od:
> od:
> for I2 from 1 to N do
> A[I2-1,N] := b[I2-1]:
> od:
> LINSYS(N,OK,A,Y):
> if OK = FALSE then
> break:
> fi:
> for I2 from 1 to N do
> K3[I2-1] := Y[I2-1]:
> X3[I2-1] := X[I2-1]+K3[I2-1]:
> od:
> for I2 from 1 to N do
> for J from 1 to N do
> A[I2-1,J-1] := evalf(P[I2,J](seq(X3[i],i=0..N-1))):
> od:
> od:
> for I2 from 1 to N do
> A[I2-1,N] := b[I2-1]:
> od:
> LINSYS(N,OK,A,Y):
> if OK = FALSE then
> break:
> fi:
> # Step 7
> for I2 from 1 to N do
> K4[I2-1] := Y[I2-1]:
> X4[I2-1] := X[I2-1]+(K1[I2-1]+2*K2[I2-1]+2*K3[I2-1]+K4[I2-1])/6:
> X[I2-1] := X4[I2-1]:
> od:
> fprintf(OUP, ` %2d`, K):
> for I2 from 1 to N do
> fprintf(OUP, ` %11.8f `, X[I2-1]):
> od:
> fprintf(OUP,` \n`):
> od:
> # Step 8
> if OUP <> default then
> fclose(OUP):
> print(`Output file `,NAME,` created sucessfully`):
> fi:
> fi: