%% ODE with Stringent Error Tolerances
% Compared to ode45, the ode113, ode78, and ode89 solvers are better at solving problems with stringent error tolerances. A common situation where these solvers excel is in orbital dynamics problems, where the solution curve is smooth and requires high accuracy in each step from the solver.
% The two-body problem considers two interacting masses m1 and m2 orbiting a common plane. In this example, one of the masses is significantly larger than the other. With the heavy body at the origin, the equations of motion are
% \begin{array}{l}
% {x^{\prime } }^{\prime } =-x/r^3 \\
% {y^{\prime } }^{\prime } =-y/r^3 ,
% \end{array}
% where
% r=\sqrt{x^2 +y^2 }\ldotp
% To solve the problem, first convert to a system of four first-order ODEs using the substitutions
% \begin{array}{l}
% y_1 =x\\
% y_2 =x^{\prime } \\
% y_3 =y\\
% y_4 =y^{\prime } \ldotp
% \end{array}
% The substitutions produce the first-order system
% \begin{array}{l}
% {y_1 }^{\prime } =y_2 \\
% {y_2 }^{\prime } =-y_1 /r^3 \\
% {y_3 }^{\prime } =y_4 \\
% {y_4 }^{\prime } =-y_3 /r^3 \ldotp
% \end{array}
% twobodyode, a local function included at the end of this example, codes the system of equations for the two-body problem.
% function dy = twobodyode(t,y)
% % Two-body problem with one mass much larger than the other.
% r = sqrt(y(1)^2 + y(3)^2);
% dy = [y(2);
%      -y(1)/r^3;
%       y(4);
%      -y(3)/r^3];
% end
% Solve the system of ODEs using ode89. Specify stringent error tolerances of 1e-13 for RelTol and 1e-14 for AbsTol.

opts = odeset('Reltol',1e-13,'AbsTol',1e-14,'Stats','on');
tspan = [0 10*pi];
y0 = [2 0 0 0.5];
[t,y] = ode89(@twobodyode, tspan, y0, opts);
plot(t,y)
legend('x','x''','y','y''','Location','SouthEast')
title('Position and Velocity Components')

%

figure
plot(y(:,1),y(:,3),'-o',0,0,'ro')
axis equal
title('Orbit of Smaller Mass')

%
% Compared to ode45, the ode89 solver is able to obtain the solution faster and with fewer steps and function evaluations.

function dy = twobodyode(t,y)
% Two-body problem with one mass much larger than the other.
r = sqrt(y(1)^2 + y(3)^2);
dy = [y(2); 
    -y(1)/r^3;
    y(4);
    -y(3)/r^3];
end

% Copyright 2021 The MathWorks, Inc.
%