(* ::Package:: *)

(************************************************************************)
(* This file was generated automatically by the Mathematica front end.  *)
(* It contains Initialization cells from a Notebook file, which         *)
(* typically will have the same name as this file except ending in      *)
(* ".nb" instead of ".m".                                               *)
(*                                                                      *)
(* This file is intended to be loaded into the Mathematica kernel using *)
(* the package loading commands Get or Needs.  Doing so is equivalent   *)
(* to using the Evaluate Initialization Cells menu command in the front *)
(* end.                                                                 *)
(*                                                                      *)
(* DO NOT EDIT THIS FILE.  This entire file is regenerated              *)
(* automatically each time the parent Notebook file is saved in the     *)
(* Mathematica front end.  Any changes you make to this file will be    *)
(* overwritten.                                                         *)
(************************************************************************)



ODE23Step[y_,f_,h_]:= Module[
{k1,k2,k3,k4,yNew,zNew},
k1=f[y];
k2=f[y+ 1/2 h k1];
k3=f[y+3/4 h k2];
yNew=y + h ( 2/9 k1 + 1/3 k2 + 4/9 k3);
k4=f[yNew];
zNew=y +h*(7/24 k1 + 1/4 k2 + 1/3 k3 + 1/8 k4);
{yNew, Norm[zNew-yNew]}
]  

DriverFixed[y0_,f_,n_, TMax_,Stepper_]:=Module[
{Data,y=y0,err,h=TMax/n,t=0.0},
Data = Table[
{y, err} = Stepper[y, f, h];
t = t + h;
{t, y},
{n}
];
(* fixes it so that it comes out in the form {ts, ys} *)
Transpose[Data]
]

DriverVariableStep1[y0_,h0_,f_,n_, ErrorOrder_,LocalTol_,Stepper_]:=Module[
{y=y0,t=0.0,h=h0,err,Data, SafetyFactor=0.9, yNew},
  Data = Table[
{yNew, err} = Stepper[y, f, h];
If[err<LocalTol, {t,y}={t+h,yNew}];
h*=SafetyFactor*(LocalTol/err)^(1/ErrorOrder);
{t, y},
{n}];
(* fixes it so that it comes out in the form {ts, ys} *)
Transpose[Data]
]


f1[{y1_,y2_,y3_,y4_,y5_,y6_}] :=Module[
{\[Mu]=1/81.45,d1,d2 },
d1 =(y2^2+y3^2+(y1+\[Mu]-1)^2)^(3/2); 
d2=(y2^2+y3^2+(y1+\[Mu])^2)^(3/2);
{y4,y5,y6,
2y5 + y1 - (\[Mu] (y1+\[Mu]-1))/d1 -((1-\[Mu])(y1+\[Mu]))/d2,
-2y4 + y2 - (\[Mu] y2)/d1 -((1-\[Mu])y2)/d2,
- ((\[Mu] y3)/d1) -((1-\[Mu])y3)/d2
}
]
IC = {0.994,0,0,0,-2.001585106370825224,0};



