(******************************************************************
 * Ordinary Differential Equation solving techniques              *
 ******************************************************************)


Ft1[t_,y_] := Return[N[y+t*(t-2.0)]];

Runge4[f_,y_,t_,h_]:=
  Block[{k1,k2,k3,k4,h2 = h/2.0},
	k1 = h*f[t,y];
	k2 = h*f[t+h2,y+k1/2.0];
	k3 = h*f[t+h2,y+k2/2.0];
	k4 = h*f[t+h ,y+k3];
	Return[y + (k1+2.0*k2+2.0*k3+k4)/6.0];];



(* y = {yn,fn,fn_1,fn_2,fn_3} *)
AdamsBashford4[f_,y0_,t_,h_]:=
  Block[{y = y0},
	y[[1]] = y[[1]] + (55.0*y[[2]]-59*y[[3]]+37*y[[4]]-9*y[[5]])*h/24.0;
	y[[5]] = y[[4]];  y[[4]] = y[[3]];  y[[3]] = y[[2]];
	y[[2]] = f[t+h,y[[1]]];
	Return[y];];

Runode[f_,y0_,t0_,tf_,h_] :=
  Block[{y=y0,t=t0},
	While[t < tf,
	      y = Runge4[f,y,t,h];
	      t = t+h;
	      Print["t = ",t," y = ",N[y,10], " ytrue = ", N[Exp[t]-t*t,10]];
	      ];];

Runode2[f_,y0_,t0_,tf_,h_] :=
  Block[{y=y0,t=t0},
	While[t < tf,
	      y = AdamsBashford4[f,y,t,h];
	      t = t+h;
	      Print["t = ",t," y = ",N[y[[1]],10], " ytrue = ", N[Exp[t]-t*t,10]];
	      ];];

(* Question 4:
   Runode[Ft1,1.0,0.0,1.0,0.1] 

t = 0.1 y = 1.095171042 ytrue = 1.095170918
t = 0.2 y = 1.181403009 ytrue = 1.181402758
t = 0.3 y = 1.25985919  ytrue = 1.259858808
t = 0.4 y = 1.331825214 ytrue = 1.331824698
t = 0.5 y = 1.398721924 ytrue = 1.398721271
t = 0.6 y = 1.462119591 ytrue = 1.4621188
t = 0.7 y = 1.523753635 ytrue = 1.523752707
t = 0.8 y = 1.585541991 ytrue = 1.585540928
t = 0.9 y = 1.649604305 ytrue = 1.649603111
t = 1.0 y = 1.718283148 ytrue = 1.718281828

Error = -1.32 x 10^(-6)
*)
   

(* Question 5:
   Runode2[Ft1,{1.0,Ft1[0,1],Ft1[-0.1,0.89483742],Ft1[-0.2,0.77873075],
		    Ft1[-0.3,0.65081822]},0.0,1.0,0.1] 

t = 0.1 y = 1.095167775 ytrue = 1.095170918
t = 0.2 y = 1.181395423 ytrue = 1.181402758
t = 0.3 y = 1.259846726 ytrue = 1.259858808
t = 0.4 y = 1.331806924 ytrue = 1.331824698
t = 0.5 y = 1.398696694 ytrue = 1.398721271
t = 0.6 y = 1.462086192 ytrue = 1.4621188
t = 0.7 y = 1.523710656 ytrue = 1.523752707
t = 0.8 y = 1.585487806 ytrue = 1.585540928
t = 0.9 y = 1.649537054 ytrue = 1.649603111
t = 1.0 y = 1.718200704 ytrue = 1.718281828

Error = 8.1124 x 10^(-5)

*)
