Not able to run

y0 = [r0 Vn0];

t0 = 0;

nrevs = 100;

tmax = nrevs*Period;

h0 = 0;

x0 = [r0; h0; V0];

r= @(x) [sqrt(x(1)^2 + x(2)^2 + x(3)^2)];

f = @(t,x) [x(3); x(2)*x(3)^2/x(1) -mu/x(1)^2; -x(2)*x(3)/x(1)];

settings1=odeset('RelTol',1e-4,'AbsTol',1e-4);

settings2=odeset('RelTol',1e-8,'AbsTol',1e-8);

settings3=odeset('RelTol',1e-11,'AbsTol',1e-11);

[T1,x1]=ode45(f,[0 tmax],x0,settings1)

[T2,x2]=ode45(f,[0 tmax],x0,settings2)

[T3,x3]=ode45(f,[0 tmax],x0,settings3)

% Plot the trajectories for each tolerance

% Plot the trajectories for each tolerance

figure;

plot(x1(:,1).*cos(x1(:,2)), x1(:,1).*sin(y1(:,2)), '-', ...

x2(:,1).*cos(x2(:,2)), x2(:,1).*sin(x2(:,2)), '-', ...

x3(:,1).*cos(x3(:,2)), x3(:,1).*sin(x3(:,2)), '-');

legend('tol=1e-4', 'tol=1e-8', 'tol=1e-11');

title('Orbit Trajectories');

Answer 1

[Sulaymon Eshkabilov]

Here is the corrected code (Note that the Initial Conditions were missing and some arbitrary values are taken. Therefore they need to be changed. Also, note that in figure 1 three solution sets overlap each other. Hence, they are plotted in separate) :

t0 = 0;
nrevs = 100;
Period = 2; 
tmax = nrevs*Period;
mu = 2.5;
h0 = 0;   % 
r0 = 1;   % use the given value
V0 = 0.5;   % use the given value
x0 = [r0; h0; V0];
f = @(t,x) ([x(3); x(2).*x(3).^2./x(1)-mu./x(1).^2; -x(2).*x(3)./x(1)]);
settings1=odeset('RelTol',1e-4,'AbsTol',1e-4);
settings2=odeset('RelTol',1e-8,'AbsTol',1e-8);
settings3=odeset('RelTol',1e-11,'AbsTol',1e-11);
 [T1,y1]=ode45(f,[0 tmax],x0,settings1);

Warning: Failure at t=2.004057e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (7.105427e-15) at time t.

 [T2,y2]=ode45(f,[0 tmax],x0,settings2);

Warning: Failure at t=2.003988e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (7.105427e-15) at time t.

 [T3,y3]=ode45(f,[0 tmax],x0,settings3);

Warning: Failure at t=2.003988e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (7.105427e-15) at time t.

 % Plot the trajectories for each tolerance
figure;
plot(y1(:,1).*cos(y1(:,2)), y1(:,1).*sin(y1(:,2)), 'r-', ...
     y2(:,1).*cos(y2(:,2)), y2(:,1).*sin(y2(:,2)), 'b--', ...
     y3(:,1).*cos(y3(:,2)), y3(:,1).*sin(y3(:,2)), 'g-.');
legend('tol=1e-4', 'tol=1e-8', 'tol=1e-11');
title('Orbit Trajectories');
axis equal  

figure
plot(y1(:,1).*cos(y1(:,2)), y1(:,1).*sin(y1(:,2)), 'r-')  
legend('tol=1e-4');
title('Orbit Trajectories');
axis equal

figure;
plot(y2(:,1).*cos(y2(:,2)), y2(:,1).*sin(y2(:,2)), 'b--');
legend('tol=1e-8');
title('Orbit Trajectories');
axis equal

figure;
plot(y3(:,1).*cos(y3(:,2)), y3(:,1).*sin(y3(:,2)), 'g-.');
legend('tol=1e-11');
title('Orbit Trajectories');
axis equal