Some corrections in the ODE function (ODAYRK1)
Tol = 1e-12;
Tol0 = 1e-9;
tspan = [0 10];
options = odeset('RelTol', Tol, 'AbsTol', Tol0);
% Assuming IDATA contains the initial conditions for your system
IDATA = [1, 0, 0, 0, 0, 0];
[t, y] = ode45(@ODAYRK1, tspan, IDATA, options);
plot(t, y), grid on
function ODAYRK = ODAYRK1(t, y)
% global M e0 mi r0
mi = 1; % <-- define it inside the ode function
% Extracting variables from the state vector
X = y(1);
Y = y(2);
Z = y(3);
Vx = y(4);
Vy = y(5);
Vz = y(6);
% Compute some intermediate values
r = sqrt(X^2 + Y^2 + Z^2);
% Define the ODEs
dXdt = Vx;
dYdt = Vy;
dZdt = Vz;
dVxdt = mi/(r^3)*X;
dVydt = mi/(r^3)*Y; % <-- correction at this line
dVzdt = mi/(r^3)*Z;
% Assemble the derivative vector
ODAYRK = [dXdt; dYdt; dZdt; dVxdt; dVydt; dVzdt];
end
