> restart:
> # BROYDEN ALGORITHM 10.2
> #
> # 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)): tolerance TOL:
> #          maximum number of iterations N.
> #
> # OUTPUT:  Approximate solution X = (X(1),...,X(n)) or a message
> #          that the number of iterations was exceeded.
> print(`This is the Broyden 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.`):
>    fi:
> od:
> for I1 from 1 to N do
>    print(`Input the function f[i](x1,x2, ... xn) in terms of x1, x2, ... for i = `):print(I1):
>    print(`For example: x1-x2^2+1`):
>    F[I1] := scanf(`%a`)[1]:print(`function is `):print(F[I1]):
> od:
> for I1 from 1 to N do
>    for J from 1 to N do
>       PD[I1,J] := diff(F[I1],evaln(x || J)):print(`partial derivative i,j `):print(I1,J):print(PD[I1,J]):
>       PD[I1,J] := unapply(PD[I1,J],evaln(x || (1..N))):
>    od:
> od:
> for I1 from 1 to N do
>    F[I1] := unapply(F[I1],evaln(x || (1..N))):
> od:
> OK := FALSE:
> while OK = FALSE do
>    print(`Input tolerance`):
>    TOL := scanf(`%f`)[1]:print(`Tolerance = `):print(TOL):
>    if TOL > 0 then
>       OK := TRUE:
>    else
>       print(`Tolerance must be positive.`):
>    fi:
> od:
> OK := FALSE:
> while OK = FALSE do
>    print(`Input the maximum number of iterations.`):
>    NN := scanf(`%d`)[1]:print(`Maximum number of iterations = `):print(NN):
>    if NN > 0 then
>       OK := TRUE:
>    else
>       print(`Must be a positive integer.`):
>    fi:
> od:
> for I1 from 1 to N do
> print(`Input initial approximation X(I), I = `, I1):
> X[I1-1] := scanf(`%f`)[1]: print(`Input is `): print(X[I1-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(`Your response 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:
>    print(`Select amount of output`):
>    print(`1. Answer only`):
>    print(`2. All intermeditate approximations`):
>    print(`Enter 1 or 2`):
>    FLAG := scanf(`%d`)[1]: print(`Your response is `):print(FLAG):
>    fprintf(OUP, `BROYDENS METHOD FOR NONLINEAR SYSTEMS\n\n`):
>    if FLAG = 2 then
>       fprintf(OUP, `Iteration, Approximation, Error\n`):
>    fi:
> # Step 1
> # A will hold the Jacobian for the initial approximation
>    for I1 from 1 to N do
>       for J from 1 to N do
>          A[I1-1,J-1] := evalf(PD[I1,J](seq(X[i],i=0..N-1))):
>       od:
> # Compute V = F(X(0))
>       V[I1-1] := evalf(F[I1](seq(X[i],i=0..N-1))):
>    od:
> # Step 2
>    for I1 from 1 to N do
>       for J from 1 to N do
>          B[I1-1,J-1] := 0:
>       od:
>       B[I1-1,I1-1] := 1:
>    od:
>    I1 := 1:
>    while I1 <= N and OK = TRUE do
>       I2 := I1+1:
>       I3 := I1:
>       if I1 <> N then
>          C := abs(A[I1-1,I1-1]):
>          for J from I2 to N do
>             if abs(A[J-1,I1-1]) > C then
>                I3 := J:
>                C := abs(A[J-1,I1-1]):
>             fi:
>          od:
>          if C <= 1.0e-20 then
>             OK := FALSE:
>          else
>             if I3 <> I1 then
>               for J from 1 to N do
>                   C := A[I1-1,J-1]:
>                   A[I1-1,J-1] := A[I3-1,J-1]:
>                   A[I3-1,J-1] := C:
>                   C := B[I1-1,J-1]:
>                   B[I1-1,J-1] := B[I3-1,J-1]:
>                   B[I3-1,J-1] := C:
>                od:
>             fi:
>          fi:
>       else
>          if abs(A[N-1,N-1]) <= 1.0e-20 then
>             OK := FALSE:
>          fi:
>       fi:
>       if OK = TRUE then
>          for J from 1 to N do
>             if J <> I1 then
>                C := A[J-1,I1-1]/A[I1-1,I1-1]:
>                for K from 1 to N do
>                   A[J-1,K-1] := A[J-1,K-1]-C*A[I1-1,K-1]:
>                   B[J-1,K-1] := B[J-1,K-1]-C*B[I1-1,K-1]:
>                od:
>             fi:
>          od:
>       fi:
>       I1 := I1+1:
>    od:
>    if OK = TRUE then
>       for I1 from 1 to N do
>          C := A[I1-1,I1-1]:
>          for J from 1 to N do
>             A[I1-1,J-1] := B[I1-1,J-1]/C:
>          od:
>       od:
>    else
>       print(`Jacobian has no inverse`):
>    fi:
>    if OK = TRUE then
> # Step 3
>       K := 2:
> # Note: S = S(1)
>       SN := 0:
>       for I1 from 1 to N do
>          S[I1-1] := 0:
>          for J from 1 to N do
>             S[I1-1] := S[I1-1]-A[I1-1,J-1]*V[J-1]:
>          od:
>          SN := SN+S[I1-1]^2:
>       od:
>       SN := sqrt(SN):
>       for I1 from 1 to N do
>          X[I1-1] := X[I1-1]+S[I1-1]:
>       od:
>       if FLAG = 2 then
>          fprintf(OUP,` %d`,K-1):
>          for I1 from 1 to N do
>             fprintf(OUP,` %11.8f`,X[I1-1]):
>          od:
>          fprintf(OUP,`\n %12.6e\n`,SN):
>       fi:
> # Step 4
>       while K < NN and OK = TRUE do
> # Step 5
>          for I1 from 1 to N do
>             VV := evalf(F[I1](seq(X[i],i=0..N-1))):
>             Y[I1-1] := VV-V[I1-1]:
>             V[I1-1] := VV:
>          od:
> # Note:  V = F(X(K)) and Y = Y(K)
> # Step 6
>          ZN := 0:
>          for I1 from 1 to N do
>             Z[I1-1] := 0:
>             for J from 1 to N do
>                Z[I1-1] := Z[I1-1]-A[I1-1,J-1]*Y[J-1]:
>             od:
>             ZN := ZN+Z[I1-1]*Z[I1-1]:
>          od:
>          ZN := sqrt(ZN):
> # Note:  Z = -A(K-1)^(-1)*Y(K)
> # Step 7
>          P := 0:
> # P will be S(K)^T*A(K-1)^(-1)*Y(K)
>          for I1 from 1 to N do
>             P := P-S[I1-1]*Z[I1-1]:
>          od:
> # Step 8
>          for I1 from 1 to N do
>             U[I1-1] := 0:
>             for J from 1 to N do
>                U[I1-1] := U[I1-1]+S[J-1]*A[J-1,I1-1]:
>             od:
> # Step 9
>          od:
>          for I1 from 1 to N do
>             for J from 1 to N do
>                A[I1-1,J-1] := A[I1-1,J-1]+(S[I1-1]+Z[I1-1])*U[J-1]/P:
>             od:
>          od:
> # Step 10
>          SN := 0:
>          for I1 from 1 to N do
>             S[I1-1] := 0:
>             for J from 1 to N do
>                S[I1-1] := S[I1-1]-A[I1-1,J-1]*V[J-1]:
>             od:
>             SN := SN+S[I1-1]^2:
>          od:
>          SN := sqrt(SN):
> # Note: A = A(K)^(-1) and S = -A(K)^(-1)*F(X(K))
> # Step 11
>          for I1 from 1 to N do
>             X[I1-1] := X[I1-1]+S[I1-1]:
>          od:
> # Note:  X = X(K+1)
>          KK := K+1:
>          if FLAG = 2 then
>             fprintf(OUP, ` %2d`, K):
>             for I1 from 1 to N do
>                fprintf(OUP, ` %11.8f`, X[I1-1]):
>             od:
>             fprintf(OUP, `\n%12.6e\n`, SN):
>          fi:
>          if SN <= TOL then
> # Procedure completed successfully
>             OK := FALSE:
>             fprintf(OUP, `Iteration number %d`, K):
>             fprintf(OUP, ` gives solution:\n\n`):
>             for I1 from 1 to N do
>                fprintf(OUP, ` %11.8f`, X[I1-1]):
>             od:
>             fprintf(OUP, `\n\nto within tolerance %.10e\n\n`, TOL):
>             fprintf(OUP, `Process is complete\n`):
>          else
> # Step 13
>             K := KK:
>          fi:
>       od:
>       if K >= NN then
> # Step 14
>          fprintf(OUP, `Procedure does not converge in %d iterations\n`, NN):
>       fi:
>    fi:
> fi:
> if OUP <> default then
> fclose(OUP):
> print(`Output file `,NAME,` created sucessfully`):
> fi: