# Solving Equations of motion; ND Solve stiff system suspected

Im modelling mercury, venus,earth and mars orbiting the sun for a 2 body problem.

Below are the initial position velocity conditions for the planets in question. Distance units are in AU Velocity units are in AU/day

Initial conditions

T = {88.0, 224.7, 365.2, 687.0};(*period in days*)
Ecc = {0.20563069, 0.00677323, 0.01671022, 0.09341233};(*eccentricity of orbit*)
\[Mu] = 2.9591220823*10^-4;(*standard gravitational parameter of sun*)
M = {174.79439000000002, 0., -2.48284000000001,
19.412480000000016};(*Mean Anomoly*)
\[Psi] = M + 2*Ecc*Sin[M Degree] + 5/4*Ecc^2*Sin[2*M Degree] +
13/12*Ecc^3*Sin[3*M Degree];(*True Anomoly*)
a = {0.38709893, 0.7233319899999999, 1.00000011, 1.52366231};
r = Table[a[[i]] (1 - Ecc[[i]]^2)/(1 + Ecc [[i]] Cos[\[Psi] [[i]] Degree ]), {i, 4}] (*orbital position*)
rx = Table[r[[i]] Cos[\[Psi][[i]] Degree], {i, 4}](*x component of position*)
ry = Table[r[[i]] Sin[\[Psi][[i]] Degree ], {i, 4}](*y component of position*)
v = Table[Sqrt[\[Mu] (2/r[[i]] - 1/a[[i]])], {i, 4}](*orbital velocity*)
vx = v Cos[\[Psi] Degree](*x component of velocity*)
vy = v Sin[\[Psi] Degree](* y component of velocity*)


Solving using NDSolve

Here i am using NDsolve to solve the ODE's for the equations of motion in x,y direction

planet = Table[NDSolve[{(x''[t] + (\[Mu]*x[t])/(x[t]^2 + y[t]^2)^(
3/2) == 0, (y''[t] + (\[Mu]*
y[t])/(x[t]^2 + y[t]^2)^(3/2) == 0, x == rx[[i]],
y == ry[[i]], x' == vx[[i]],
y' == vy[[i]]}, {x, y}, {t,800, 0.1}], {i,4}];


The problem

If i were to set the intial conditons to the magnitudes of postions and velocity, the simulation works fine, however, its when i use the components of (x,y) for position and velocity that i encounter the error message.

NDSolve::ndsz: At t == 76.8007476284912, step size is effectively zero; singularity or stiff system suspected.

NDSolve::ndsz: At t == 204.5293849935109, step size is effectively zero; singularity or stiff system suspected.

NDSolve::ndsz: At t == 333.0376141002715, step size is effectively zero; singularity or stiff system suspected.

General::stop: Further output of NDSolve::ndsz will be suppressed during this calculation.

What can i do to fix this?

It is necessary to correctly calculate the orbital velocity

T = {88.0, 224.7, 365.2, 687.0};(*period in days*)Ecc = {0.20563069,
0.00677323, 0.01671022,
0.09341233};(*eccentricity of orbit*)\[Mu] = \
0.0002959122082834964543838926328576059910.122911634152604;\
(*standard gravitational parameter of sun*)M = {174.79439000000002,
0., -2.48284000000001,
19.412480000000016};(*Mean Anomoly*)\[Psi] =
M + 2*Ecc*Sin[M Degree] + 5/4*Ecc^2*Sin[2*M Degree] +
13/12*Ecc^3*Sin[3*M Degree];(*True Anomoly*)a = {0.38709893,
0.7233319899999999, 1.00000011, 1.52366231};
r = Table[
a[[i]] (1 - Ecc[[i]]^2)/(1 + Ecc[[i]] Cos[\[Psi][[i]] Degree]), {i,
4}] (*orbital position*);
rx = Table[
r[[i]] Cos[\[Psi][[i]] Degree], {i, 4}](*x component of position*);
ry = Table[
r[[i]] Sin[\[Psi][[i]] Degree], {i, 4}](*y component of position*);
v = Table[
Sqrt[\[Mu] (2/r[[i]] - 1/a[[i]])], {i, 4}](*orbital velocity*);
vx = -v Sin[\[Psi] Degree](*x component of velocity*);
vy = v Cos[\[Psi] Degree](*y component of velocity*);

pl = {Mercury, Venus, Earth, Mars};

planet = Table[
NDSolve[{x''[t] + (\[Mu]*x[t])/(x[t]^2 + y[t]^2)^(3/2) == 0,
y''[t] + (\[Mu]*y[t])/(x[t]^2 + y[t]^2)^(3/2) == 0,
x == rx[[i]], y == ry[[i]], x' == vx[[i]],
y' == vy[[i]]}, {x, y}, {t, 0, 800}], {i, 4}];

{Table[ParametricPlot[{x[t], y[t]} /. planet[[i]], {t, 0, 800},
PlotLabel -> pl[[i]]], {i, 1, Length[planet]}],
ParametricPlot[
Evaluate[Table[{x[t], y[t]} /. planet[[i]], {i, Length[pl]}]], {t,
0, 800}, PlotLegends -> pl]} 