Gravitational Problem: p29

f1[{y1_,y2_,y3_,y4_,y5_,y6_},μ_] :={y4,y5,y6,
2y5 + y1 - ((μ(y1+μ-1))/((((y2^2)+(y3^2)+(((y1+μ-1))^2)))^(3/2))) -(((1-μ)(y1+μ))/((((y2^2)+(y3^2)+(((y1+μ))^2)))^(3/2))),
-2y4 + y2 - ((μ y2)/((((y2^2)+(y3^2)+(((y1+μ-1))^2)))^(3/2))) -(((1-μ)y2)/((((y2^2)+(y3^2)+(((y1+μ))^2)))^(3/2))),
- ((μ y3)/((((y2^2)+(y3^2)+(((y1+μ-1))^2)))^(3/2))) -(((1-μ)y3)/((((y2^2)+(y3^2)+(((y1+μ))^2)))^(3/2)))
}
IC = {0.994,0,0,0,-2.001585106370825224,0};
TMax=17.06521656015796; 
μVal =1/81.45;
ySol[t_]=y[t]/.NDSolve[{
y'[t]==f1[y[t],μVal],
y[0]==IC
},
y[t],
{t, 0, TMax}
][[1]];
Plot[ySol[t],{t, 0, TMax}]
ParametricPlot[ Take[ ySol[t],2],{t,0,TMax}]

TableForm[Transpose[{IC, ySol[TMax]}]]
IC-ySol[TMax]

Gravitational Problem: p30

f1[{y1_,y2_,y3_,y4_,y5_,y6_},μ_] :={y4,y5,y6,
2y5 + y1 - ((μ(y1+μ-1))/((((y2^2)+(y3^2)+(((y1+μ-1))^2)))^(3/2))) -(((1-μ)(y1+μ))/((((y2^2)+(y3^2)+(((y1+μ))^2)))^(3/2))),
-2y4 + y2 - ((μ y2)/((((y2^2)+(y3^2)+(((y1+μ-1))^2)))^(3/2))) -(((1-μ)y2)/((((y2^2)+(y3^2)+(((y1+μ))^2)))^(3/2))),
- ((μ y3)/((((y2^2)+(y3^2)+(((y1+μ-1))^2)))^(3/2))) -(((1-μ)y3)/((((y2^2)+(y3^2)+(((y1+μ))^2)))^(3/2)))
}
IC = {0.87978,0,0,0,-0.3797,0};
TMax=19.140457061662071; 
μVal =1/81.45;
ySol[t_]=y[t]/.NDSolve[{
y'[t]==f1[y[t],μVal],
y[0]==IC
},
y[t],
{t, 0, TMax}
][[1]];
Plot[ySol[t],{t, 0, TMax}]
ParametricPlot[ Take[ ySol[t],2],{t,0,TMax}]

TableForm[Transpose[{IC, ySol[TMax]}]]
IC-ySol[TMax]

Euler

f[y_]:=(y^2);
h=0.1;
{tt,yy}={0.0,1.0};
Data=Table[
{tt,yy}={tt+h, yy+ h* f[yy]},
{1/h}
]
Show[
Plot[(1/(1-t)),{t, 0, 0.99},
AxesLabel->{"t","y"}],
ListPlot[Data]
]

f[y_]:=(y^2);
EulerVal[n_]:= Module[
{tt=0.0,yy=1.0,h=0.5/n},
Data=Table[
{tt,yy}={tt+h, yy+ h* f[yy]},
{n}
][[-1]]
]
EulerVal[120]

ListPlot[
Table[{(0.5/n),EulerVal[n][[2]]},{n, 1, 542}],
AxesLabel->{"h","y"}
]

ListPlot[
Table[{(0.5/n),2-EulerVal[n][[2]]},{n, 1, 542}],
AxesLabel->{"h","Absolute Error t=0.5"}
]