(*
*  MULLER'S ALGORITHM 2.8
*
*  To find a solution to the equation f(x) = 0 
*  given three approximations x0, x1 and x2:
*
*  INPUT: x0, x1, x2; tolerance TOL;
*      maximum number of iterations NO.
*
*  OUTPUT: approximate solution p or a message that the
*      method fails.
*
*  This implementation allows for a switch to complex arithmetic.
*  The coefficients are stored in the vector A, so the dimension
*  of A may have to be changed.
*)
TEMP = Input["This is Muller's Method.\n
Input the function P(X) in terms of x.\n\n
For example: x^2 - 2x + 2 \n"];
P[x_] := Evaluate[TEMP];
OK = 1;
If[Exponent[P[x],x] == 0,
   Print["Polynomial is constant function."];
   OK = 0;
];
If[Exponent[P[x],x] == 1,
   P = -Coefficient[P[x],x,0]/Coefficient[P[x],x,1];
   Print["Polynomial is linear: root is ",SetPrecision[P,10]];
   OK = 0;
];
If[OK == 1,
   OK = 0;    
   While[OK == 0,
      TOL=Input["Input the tolerance\n"];
      If[TOL<=0,
         Input["Tolerance must be positive.\n
         Enter 1 [enter] to continue\n"],
         OK = 1;
      ];
   ];
   OK = 0;
   While[OK == 0,
      M = Input["Input maximum number of iterations\n
        no decimal point\n"];
      If[M<=0,
         Input["Must be a positive integer.\n
         Enter 1 [enter] to continue\n"],
         OK = 1;
      ];
   ];
   X[0]=Input["Input the first of three starting values\n"]; 
   X[1]=Input["Input the second of three starting values\n"];
   X[2]=Input["Input the third starting value\n"];
];
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,"Muller's Method\n"];
      Write[OUP,"The output is i, approximation x(i), f(x(i))\n\n"];
      Write[OUP,"Three columns means the results are real numbers,\n"];
      Write[OUP,"Five columns means the results are complex\n"];
      Write[OUP,"numbers with real and imaginary parts of x(i)\n"];
      Write[OUP,"followed by real and imaginary parts of f(x(i)).\n"];
      Write[OUP,"\n"];
   ];
   F[0]=P[X[0]];
   F[1]=P[X[1]];
   F[2]=P[X[2]];
(*STEP 1*)
   H[0]=X[1]-X[0];
   H[1]=X[2]-X[1];
   DEL1[0]=(F[1]-F[0])/H[0];
   DEL1[1]=(F[2]-F[1])/H[1];
   DEL=(DEL1[1]-DEL1[0])/(H[1]+H[0]);
   i=3;
(*STEP 2*)
   While[i<= M && OK == 1,
(*STEP 3*)
      B=DEL1[1]+H[1]*DEL;
      d=B*B-4*F[2]*DEL;
(* test to see if straight line *)
      If[Abs[DEL]<= 10^-20,
(* straight line - test if horizontal line*)
         If[Abs[DEL1[1]]<= 10^-20,
            Print["Horizontal Line"];
            OK = 0,
(* straight line but not horizontal *)
            X[3]=(F[2]-DEL1[1]*X[2])/DEL1[1];
            H[2]=X[3]-X[2];
         ],
(* not a straight line *)
         d=Sqrt[d];
(*STEP 4*)
         e=B+d;
         If[Abs[B-d]>Abs[e],
            e=B-d;
         ];
(*STEP 5*)
         H[2]=-2*F[2]/e;
         X[3]=X[2]+H[2];
      ];
      If[OK == 1,
         F[3] = N[P[X[3]]];
    Write[OUP,i," ",SetPrecision[X[3],10],"   ",SetPrecision[F[3],10]];
      ];
(*STEP 6*)
      If[Abs[H[2]]<TOL,
         Write[OUP,"Method succeeds\n"];
         Write[OUP,"Approximation is within \n",SetPrecision[TOL,10]];
         Write[OUP,"Number of iterations = \n",i];   
         OK=0,
(*STEP 7*)
         For[J = 1,
            J <= 2,
            J++,
            H[J-1]=H[J];
            X[J-1]=X[J];
            F[J-1]=F[J];
         ];
         X[2]=X[3];
         F[2]=F[3];
         DEL1[0]=DEL1[1];
         DEL1[1]=(F[2]-F[1])/H[1];
         DEL=(DEL1[1]-DEL1[0])/(H[1]+H[0]);
      ];
      i=i+1;
   ];
(*STEP 8*)
   If[i > M && OK ==1,
      Print["Method failed\n"];
   ];
   If[OUP == "OutputStream[",NAME," 3]",
      Print["Output file: ",NAME," created successfully\n"];
      Close[OUP];
   ];
];
Quit[];