> 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)). 
> #          
> alg104 := proc() local LINSYS,OK,N,I,F,J,P,NN,X,FLAG,NAME,OUP,K,A,b,X1,X2,X3,X4,K1,K2,K3,K4,H;
> LINSYS := proc(N,OK,A,Y) local K, I, Z, IR, IA, J, C, L, JA;
> K := N-1;
> OK := TRUE;
> I := 1;
> while OK = TRUE and I <= K do
> Z := abs(A[I-1,I-1]);
> IR := I;
> IA := I+1;
> for J from IA to N do
> if abs(A[J-1,I-1]) > Z then
> IR := J;
> Z := abs(A[J-1,I-1]);
> fi;
> od;
> if Z <= 1.0e-20 then
> OK := FALSE;
> else
> if IR <> I then
> for J from I to N+1 do
> C := A[I-1,J-1];
> A[I-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,I-1]/A[I-1,I-1];
> if abs(C) <= 1.0e-20 then
> C := 0;
> fi;
> for L from I to N+1 do
> A[J-1,L-1] := A[J-1,L-1]-C*A[I-1,L-1];
> od;
> od;
> fi;
> I := I+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 I from 1 to K do
> J := N-I;
> 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
> printf(`Linear system is singular\n`);
> fi;
> end;
> printf(`This is the Continuation Method for Nonlinear Systems.\n`);
> OK := FALSE;
> while OK = FALSE do
> printf(`Input the number n of equations.\n`);
> N := scanf(`%d`)[1];
> if N >= 2 then
> OK := TRUE;
> else
> printf(`N must be an integer greater than 1.\n`);
> fi;
> od;
> for I from 1 to N do
> printf(`Input the function F%d in terms of x1 ... x%d\n` ,I,N);
> F[I] := scanf(`%a`)[1];
> od;
> for I from 1 to N do
> for J from 1 to N do
> P[I,J] := diff(F[I],evaln(x . J));
> P[I,J] := unapply(P[I,J],evaln(x . (1..N)));
> od;
> od;
> for I from 1 to N do
> F[I] := unapply(F[I],evaln(x . (1..N)));
> od;
> OK := FALSE;
> while OK = FALSE do
> printf(`Input the number N for RK4.\n`);
> NN := scanf(`%d`)[1];
> if NN > 0 then
> OK := TRUE;
> else
> printf(`Must be a positive integer.\n`);
> fi;
> od;
> for I from 1 to N do
> printf(`Input initial approximation X(%d).\n`, I);
> X[I-1] := scanf(`%f`)[1];
> od;
> if OK = TRUE then
> printf(`Select output destination\n`);
> printf(`1. Screen\n`);
> printf(`2. Text file\n`);
> printf(`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;
> # Step 1
> H := 1/NN;
> for I from 1 to N do
> b[I-1] := H*evalf(-F[I](seq(X[i],i=0..N-1)));
> od;
> # Step 2
> for K from 1 to NN do
> # Steps 3 - 6
> for I from 1 to N do
> for J from 1 to N do
> A[I-1,J-1] := evalf(P[I,J](seq(X[i],i=0..N-1)));
> od;
> od;
> for I from 1 to N do
> A[I-1,N] := b[I-1];
> od;
> LINSYS(N,OK,A,Y);
> if OK = FALSE then
> break;
> fi;
> for I from 1 to N do
> K1[I-1] := Y[I-1];
> X1[I-1] := X[I-1]+0.5*K1[I-1];
> od;
> for I from 1 to N do
> for J from 1 to N do
> A[I-1,J-1] := evalf(P[I,J](seq(X1[i],i=0..N-1)));
> od;
> od;
> for I from 1 to N do
> A[I-1,N] := b[I-1];
> od;
> LINSYS(N,OK,A,Y);
> if OK = FALSE then
> break;
> fi;
> for I from 1 to N do
> K2[I-1] := Y[I-1];
> X2[I-1] := X[I-1]+0.5*K2[I-1];
> od;
> for I from 1 to N do
> for J from 1 to N do
> A[I-1,J-1] := evalf(P[I,J](seq(X2[i],i=0..N-1)));
> od;
> od;
> for I from 1 to N do
> A[I-1,N] := b[I-1];
> od;
> LINSYS(N,OK,A,Y);
> if OK = FALSE then
> break;
> fi;
> for I from 1 to N do
> K3[I-1] := Y[I-1];
> X3[I-1] := X[I-1]+K3[I-1];
> od;
> for I from 1 to N do
> for J from 1 to N do
> A[I-1,J-1] := evalf(P[I,J](seq(X3[i],i=0..N-1)));
> od;
> od;
> for I from 1 to N do
> A[I-1,N] := b[I-1];
> od;
> LINSYS(N,OK,A,Y);
> if OK = FALSE then
> break;
> fi;
> # Step 7
> for I from 1 to N do
> K4[I-1] := Y[I-1];
> X4[I-1] := X[I-1]+(K1[I-1]+2*K2[I-1]+2*K3[I-1]+K4[I-1])/6;
> X[I-1] := X4[I-1];
> od;
> fprintf(OUP, ` %2d`, K);
> for I from 1 to N do
> fprintf(OUP, ` %11.8f `, X[I-1]);
> od;
> fprintf(OUP,` \n`);
> od;
> # Step 8
> if OUP <> default then
> fclose(OUP):
> printf(`Output file %s created sucessfully`,NAME);
> fi;
> fi;
> RETURN(0);
> end;
> alg104();