> restart:
> # NONLINEAR SHOOTING ALGORITHM 11.2
> #
> # To approximate the solution of the nonlinear boundary-value problem
> #
> #          Y'' = F(X,Y,Y'), A<=X<=B, Y(A) = ALPHA, Y(B) = BETA:
> #
> #
> # INPUT:   Endpoints A,B: boundary conditions ALPHA, BETA: number of
> #          subintervals N: tolerance TOL: maximum number of iterations M.
> #
> # OUTPUT:  Approximations W(1,I) TO Y(X(I)): W(2,I) TO Y'(X(I))
> #          for each I=0,1,...,N or a message that the maximum
> #          number of iterations was exceeded.
> print(`This is the Nonlinear Shooting Method.`):
> print(`Input the function F(X,Y,Z) in terms of x, y, z.`):
> print(`For example:   (32+2*x^3-y*z)/8`):
> F := scanf(`%a`)[1]: print(`F(x,y,z) = `):print(F):
> FY := diff(F,y):
> FYP := diff(F,z):
> F := unapply(F,x,y,z):
> FY := unapply(FY,x,y,z):
> FYP := unapply(FYP,x,y,z):
> OK := FALSE:
> while OK = FALSE do
> print(`Input left and right endpoints separated by blank.`):
> A := scanf(`%e`)[1]: print(`Left endpoint = `):print(A):
> B := scanf(`%e`)[1]: print(`Right endpoint = `):print(B):
> if A >= B then
> print(`Left endpoint must be less than right endpoint.`):
> else OK := TRUE:
> fi:
> od:
> print(`Input Y(A).`):
> ALPHA := scanf(`%e`)[1]: print(`alpha = `):print(ALPHA):
> print(`Input Y(B).`): 
> BETA := scanf(`%e`)[1]: print(`beta = `):print(BETA):
> TK := (BETA-ALPHA)/(B-A):
> print(`TK := `): print(TK):
> print(`Input new TK? Enter 1 for yes or 2 for no.`):
> AA := scanf(`%d`)[1]: print(`Input is `):print(AA):
> if AA = 1 then
> print(`Input new TK`):
> TK := scanf(`%f`)[1]: print(`New TK = `):print(TK):
> fi:
> OK := FALSE:
> while OK = FALSE do
> print(`Input an integer > 1 for the number of subintervals.`):
> N := scanf(`%d`)[1]: print(`n = `):print(N):
> if N <= 1 then
> print(`Number must exceed 1.`):
> else
> OK := TRUE:
> fi:
> od:
> OK := FALSE:
> while OK = FALSE do
> print(`Input Tolerance.`):
> TOL := scanf(`%f`)[1]: print(`Tolerance = `):print(TOL):
> if TOL <= 0 then
> print(`Tolerance must be positive.`):
> else
> OK := TRUE:
> fi:
> od:
> OK := FALSE:
> while OK = FALSE do
> print(`Input maximum number of iterations.`):
> NN := scanf(`%d`)[1]:print(`Maximum number of iterations = `):print(NN):
> if NN <= 0 then
> print(`Must be positive integer.`):
> else
> OK := TRUE:
> fi:
> od:
> if OK = TRUE then
> print(`Choice of output method:`):
> print(`1. Output to screen`):
> print(`2. Output to text File`):
> print(`Please 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, `NONLINEAR SHOOTING METHOD\n\n`):
> fprintf(OUP, `  I    X(I)         W1(I)        W2(I)\n`):
> # Step 1
> H := (B-A)/N:
> K := 1:
> # TK is already computed.
> OK := FALSE:
> # Step 2
> while K <= NN and OK = FALSE do
> # Step 3
> W1[0] := ALPHA:
> W2[0] := TK:
> U1 := 0 :
> U2 := 1:
> # Step 4
> # Runge-Kutta method for systems is used in Steps 5 and 6
> for I2 from 1 to N do
> # Step 5
> X := A+(I2-1)*H:
> T := X+0.5*H:
> # Step 6
> K11 := H*W2[I2-1]:
> K12 := H*F(X,W1[I2-1],W2[I2-1]):
> K21 := H*(W2[I2-1]+0.5*K12):
> K22 := H*F(T,W1[I2-1]+0.5*K11,W2[I2-1]+0.5*K12):
> K31 := H*(W2[I2-1]+0.5*K22):
> K32 := H*F(T,W1[I2-1]+0.5*K21,W2[I2-1]+0.5*K22):
> K41 := H*(W2[I2-1]+K32):
> K42 := H*F(X+H,W1[I2-1]+K31,W2[I2-1]+K32):
> W1[I2] := W1[I2-1]+(K11+2*(K21+K31)+K41)/6:
> W2[I2] := W2[I2-1]+(K12+2*(K22+K32)+K42)/6:
> K11 := H*U2:
> K12 := H*(FY(X,W1[I2-1],W2[I2-1])*U1+FYP(X,W1[I2-1],W2[I2-1])*U2):
> K21 := H*(U2+0.5*K12):
> K22 := H*(FY(T,W1[I2-1],W2[I2-1])*(U1+0.5*K11)+FYP(T,W1[I2-1],W2[I2-1])*(U2+0.5*K21)):
> K31 := H*(U2+0.5*K22):
> K32 := H*(FY(T,W1[I2-1],W2[I2-1])*(U1+0.5*K21)+FYP(T,W1[I2-1],W2[I2-1])*(U2+0.5*K22)):
> K41 := H*(U2+K32):
> K42 := H*(FY(X+H,W1[I2-1],W2[I2-1])*(U1+K31)+FYP(X+H,W1[I2-1],W2[I2-1])*(U2+K32)):
> U1 := U1+(K11+2*(K21+K31)+K41)/6:
> U2 := U2+(K12+2*(K22+K32)+K42)/6:
> od:
> # Step 7
> # Test for accuracy
> if abs(W1[N]-BETA) < TOL then
> # Step 8
> I2 := 0:
> fprintf(OUP, `%3d %13.8f %13.8f %13.8f\n`, I2, A, ALPHA, TK):
> for I2 from 1 to N do
> J := I2+1:
> X := A+I2*H:
> fprintf(OUP, `%3d %13.8f %13.8f %13.8f\n`, I2, X, W1[J-1], W2[J-1]):
> od:
> fprintf(OUP, `Convergence in %d iterations\n`, K):
> fprintf(OUP, ` t = %14.7e\n`, TK):
> # Step 9
> OK := TRUE:
> else
> # Step 10
> # Newton's method applied to improve TK
> TK := TK-(W1[N]-BETA)/U1:
> K := K+1:
> fi:
> od:
> # Step 11
> # Procedure completed unsuccessfully.
> if OK = FALSE then
> fprintf(OUP, `Method failed after %d iterations\n`, NN):
> fi:
> fi:
> if OUP <> default then
> fclose(OUP):
> print(`Output file `,NAME,` created successfully`):
> fi: