(**********************************************************************
 *
 * Assignment 3, Mullers method
 * Questions 1,2,3
 **********************************************************************)


Mullerstep[f_,x0_,x1_,x2_]:=
  Block[{x3,a,b,f11,f12,f1,f2,det,d1,d2},
	f2 = f[x2]; f1 = f[x1];
	f11 = (f[x0]-f1)/(x0 - x1);
	f12 = (f1 -  f2)  /(x1 - x2);
	a   = (f11 - f12) /(x0 - x2);
	b   = f12 + a * (x2 - x1);
	det = N[Sqrt[b*b -4*a*f2]];
	d1 = N[b+det]; d2 = N[b-det];
	N[If[Abs[d1] >= Abs[d2], Return[x2 - 2*f2/d1], Return[x2 - 2*f2/d2]]]
];

Muller[f_,xa_,eps_]:=
  Block[{x0 = xa-0.1, x1 = xa, x2 = xa+0.1, x3},
	While[Abs[f[x2]] > eps,
	      x3 = N[Mullerstep[f,x0,x1,x2]];
	      Print[N[x3,10]];
	      x0 = x1; x1 = x2; x2 = x3;
	      ];
	Return[N[x2,10]]
	];



F1[x_]:= Return[x*Sin[x] - Cos[x]];
F2[x_]:= Return[F1[x]/(x-0.86033333)];
F3[x_]:= Return[F2[x]/(x-3.42561846)];
F4[x_]:= Return[F3[x]/(x-6.43729817)];

(* the above three numbers : 0.86033333, 3.4256184, 6.43729817 are the
first 3 positive roots of x = tan x; They were found by successively
convolving the original function with the roots found to remove
convergence to them, this allows specification of initial points far
from the roots, here the initial points used were 2.0, 2.0, 3.0 for
the three roots *)

(* Bessels first function *)

J1[x_] := Block[{k = 36, i, f=0,a},
		If[x == 0,Return[0.0],];
	       a = N[-x*x/4.0];
	       For[i=0, i <= k, i++,
		   f = N[f + a^i/(i! * (i+1)!)]];
	       Return[x*0.5*f];
	       ];

J1b[x_] := Return[J1[x]/(x-3.83170597)];
J1c[x_] := Return[J1b[x]/(x-7.01558667)];
J1d[x_] := Return[J1c[x]/(x-10.17346814)];


(* Roots using Bessel1, k = 3
   initial 2.0, root -> 3.6578901

   Using k = 4,
   3.869330992
   5.169438

   Using k = 5,
   3.828731311

 k = 6
 3.831905947
 6.42946

 k = 12
 3.83170597
 7.01558355
 10.27575119

 k = 18
 3.83170597
 7.01558667
 10.17346818
  13.32185771

 k = 26
 3.83170597
 7.01558667
 10.17346814
 13.32369194

 k = 36
 3.83170597
 7.01558667
 10.17346814
 13.32369194

 Hence the solution obtained for k = 36 are converged to the no. of digits shown
 *)
