(*  TRAPEZOIDAL WITH NEWTON ITERATION ALGORITHM 5.8
*
*   To approximate the solution of the initial value problem:
*    Y' = F(T,Y), A <= T <= B, Y(A) = ALPHA,
*    at (n+1) equally spaced numbers in the interval [A,B].
*
*   INPUT: endpoints a, b; initial condition alpha; integer n;
*           tolerance TOL; maximum number fo iterations
*	    at any one step M.
*  
*  OUTPUT: approximation W to Y at the (n+1) values of T
*	    or a message of failure.
*)
TEMP = Input["This is the Implicit Trapezoidal Method\n
 Input the function F(T,Y) in terms of t and y\n 
 \n
 For example: y-t^2+1\n"];
F[t_,y_] := Evaluate[TEMP];
DER = D[TEMP,y];
FYP[t_,y_] := Evaluate[DER];
OK = 0;
While[OK == 0,
   A = Input["Input the left endpoint\n"];
   B = Input["Input the right endpoint\n"];
   If[A >= B,
      Input["Left endpoint must be less than right endpoint\n
      \n
      Press 1 [enter] to continue\n"],
      OK = 1;
   ];
];
ALPHA = Input["Input the initial condition\n"];
OK = 0;
While[OK == 0,
   n = Input["Input a positive integer for the number of\n
   subintervals\n"];
   If[n <= 0,
      Input["Number must be a positive integer\n
      \n
      Press 1 [enter] to continue\n"],
      OK = 1;
   ];
];
OK = 0;
While[OK == 0,
   TOL = Input["Input tolerance"];
   If[TOL <= 0,
      Input["Tolerance must be positive\n
      \n
      Press 1 [enter] to continue\n"],
      OK = 1;
   ];
];
OK = 0;
While[OK == 0,
   M = Input["Input maximum number of iterations"];
   If[M <= 0,
      Input["Number of iterations must be positive\n
      \n
      Press 1 [enter] to continue\n"],
      OK = 1;
   ];
];
If[OK == 1,
   FLAG = Input["Select output destination\n
                 1. Screen\n
                 2. Text file\n
                 Enter 1 or 2\n"];
   If[FLAG == 2,
      NAME = InputString["Input the file name\n
                          For example:   output.dta\n"];
      OUP = OpenWrite[NAME,FormatType->OutputForm],
      OUP = "stdout";
   ];
   Write[OUP,"IMPLICIT TRAPEZOIDAL METHOD USING NEWTONS METHOD\n"];
   Write[OUP,"\n"];
   Write[OUP,"t         w         #iter\n"];
   Write[OUP,"\n"]; 
   (* Step 1 *) 
   T = N[A];
   W = ALPHA;
   H = (B-A)/n;
   Write[OUP,SetPrecision[T,9],"       ",SetPrecision[W,9]
         ,"        0\n"];
   i = 1;
   OK = 1;
   (* Step 2 *)
   While[i <= n && OK == 1,
      (* Step 3 *)
      XK1 = W+0.5*H*N[F[T,W]];
      W0 = XK1;
      J = 1;
      DONE = 0;
      (* Step 4 *)
      While[DONE == 0 && OK == 1,
         (* Step 5 *)
	 W = W0-(W0-XK1-0.5*H*N[F[T+H,W0]])
                  /(1-0.5*H*N[FYP[T+H,W0]]);
         (* Step 6 *)
         If[Abs[W-W0] < TOL,
	    DONE = 1;
            (* Step 7 *)
	    T = N[A+i*H];
	    Write[OUP,SetPrecision[T,9],"       ",SetPrecision[W,9],
                  "       ",J];
	    i = i+1,
	    J = J+1;
	    W0 = W;
	    If[ J > M,
	       OK = 0;
	    ];
         ];
      ];
   ];
   If[OK == 0,
      Write[OUP,"Maximum number of iterations exceeded\n"];
   ];
   (* Step 8 *)
   If[OUP == "OutputStream[",NAME," 3]",
      Print["Output file: ",NAME," created successfully\n"];
      Close[OUP];
   ];
];
Quit[];