> restart:
> # LINEAR SHOOTING ALGORITHM 11.1
> #
> # To approximate the solution of the boundary-value problem
> #
> # -Y'' + P(X)Y' + Q(X)Y + R(X) = 0, A<=X<=B, Y(A)=ALPHA, Y(B)=BETA:
> #
> #
> # INPUT: Endpoints A,B: boundary conditions ALPHA, BETA: number of
> #        subintervals N.
> #
> # OUTPUT: Approximations W(1,I) to Y(X(I)): W(2,I) to Y'(X(I))
> #         for each I=0,1,...,N.
> print(`This is the Linear Shooting Method.`):
> print(`Input the functions P(X), Q(X) and R(X) in terms of x, separated by spaces.`):
> print(`For example: -2/x 2/(x^2) sin(log(x))/(x^2)`):
> P := scanf(`%a`)[1]: print(`P(x) = `):print(P):
> Q := scanf(`%a`)[1]: print(`Q(x) = `):print(P):
> R := scanf(`%a`)[1]: print(`R(x) = `):print(R):
> P := unapply(P,x):
> Q := unapply(Q,x):
> R := unapply(R,x):
> OK := FALSE:
> while OK = FALSE do
> print(`Input left and right endpoints separated by blank.`):
> A := scanf(`%f`)[1]: print(`Left endpoint = `):print(A):
> B := scanf(`%f`)[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(`%a`)[1]: print(`alpha = `):print(ALPHA):
> print(`Input Y(B).`): 
> BETA := scanf(`%a`)[1]: print(`beta = `):print(BETA):
> OK := FALSE:
> while OK = FALSE do
> print(`Input a positive integer for the number of subintervals.`):
> N := scanf(`%d`)[1]: print(`Input = `): print(N):
> if N <= 0 then
> print(`Number must be a 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 = `): 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, `LINEAR SHOOTING METHOD\n\n`):
> fprintf(OUP, `  I  X(I)        W(1,I)      W(2,I)\n`):
> # Step 1
> H := (B-A)/N:
> U1 := ALPHA:
> U2 := 0:
> V1 := 0:
> V2 := 1:
> # Step 2
> for I2 from 1 to N do
> # Step 3
> X := A+(I2-1)*H:
> T := X+0.5*H:
> # Step 4
> K11 := H*U2:
> K12 := H*(P(X)*U2+Q(X)*U1+R(X)):
> K21 := H*(U2+0.5*K12):
> K22 := H*(P(T)*(U2+0.5*K12)+Q(T)*(U1+0.5*K11)+R(T)):
> K31 := H*(U2+0.5*K22):
> K32 := H*(P(T)*(U2+0.5*K22)+Q(T)*(U1+0.5*K21)+R(T)):
> T := X+H:
> K41 := H*(U2+K32):
> K42 := H*(P(T)*(U2+K32)+Q(T)*(U1+K31)+R(T)):
> U1 := U1+(K11+2*(K21+K31)+K41)/6:
> U2 := U2+(K12+2*(K22+K32)+K42)/6:
> K11 := H*V2:
> K12 := H*(P(X)*V2+Q(X)*V1):
> T := X+0.5*H:
> K21 := H*(V2+0.5*K12):
> K22 := H*(P(T)*(V2+0.5*K12)+Q(T)*(V1+0.5*K11)):
> K31 := H*(V2+0.5*K22):
> K32 := H*(P(T)*(V2+0.5*K22)+Q(T)*(V1+0.5*K21)):
> T := X+H:
> K41 := H*(V2+K32):
> K42 := H*(P(T)*(V2+K32)+Q(T)*(V1+K31)):
> V1 := V1+(K11+2*(K21+K31)+K41)/6:
> V2 := V2+(K12+2*(K22+K32)+K42)/6:
> U[0,I2-1] := U1:
> U[1,I2-1] := U2:
> V[0,I2-1] := V1:
> V[1,I2-1] := V2:
> od:
> # Step 5
> W1 := ALPHA:
> Z := (BETA-U[0,N-1])/V[0,N-1]:
> X := A:
> I2 := 0:
> fprintf(OUP, `%3d %11.8f %11.8f %11.8f\n`, I2, X, W1, Z):
> # Step 6
> for I2 from 1 to N do
> X := A+I2*H:
> W1 := U[0,I2-1]+Z*V[0,I2-1]:
> W2 := U[1,I2-1]+Z*V[1,I2-1]:
> fprintf(OUP, `%3d %11.8f %11.8f %11.8f\n`, I2, X, W1, W2):
> od:
> # Step 7
> if OUP <> default then
> fclose(OUP):
> print(`Output file `,NAME,` created successfully`):
> fi:
> fi: