clc
clear all
format long

IC = [1,0];
TMax = 2;
n = 10000;

h = (TMax/(n));
tolerance  = [1e-04,1e-06,1e-08,1e-10];

for i = 1:1:4
    [tt,yy,steps,res] = Driver(IC,n,h,TMax,tolerance(i));
    steps=steps+5;
    average_iteration = mean(steps);
    iter(i) = average_iteration;
    R(i) = mean(res);
end

yy1 = yy(1:1:n,1);
yy2 = yy(1:1:n,2);

h1 = figure(1)
plot(yy1,yy2)
title('Gravitational problem solution')
xlabel('y1')
ylabel('y2')
grid on
saveas(h1,'solution.fig')

h2 = figure(2)
semilogx(tolerance,iter,'-r*')
title('Evaluation')
xlabel('Tolerance')
ylabel('Average Iterations')
grid on
saveas(h2,'iterations.fig')

h3 = figure(3)
semilogx(tolerance,R,'-r*')
title('Evaluation')
xlabel('Tolerance')
ylabel('Residual R')
grid on
saveas(h3,'residual.fig')

h4 = figure(4)
plot(0:h:tt,res,'-r')
title('Evaluation')
xlabel('Time')
ylabel('Residual R')
grid on
saveas(h4,'residual_time.fig')

% %% Error Order Check
% 
% h = [((2*pi)/(n)),((2*pi)/(2*n)),((2*pi)/(4*n)),((2*pi)/(8*n)),((2*pi)/(16*n))];
%  
% %% Analytical Solution
% options = odeset('RelTol',1e-12,'AbsTol',[1e-12 1e-12 1e-12 1e-12]);
% [T,Y] = ode45(@orbit,[0,TMax],[1-e,0,0,sqrt((1+e)/(1-e))],'options');
% size_ode = size(Y,1);
% 
% %%
% for i  = 1:1:5
%     [tt,yy,steps] = Driver(IC,n,h(i),TMax);
%     size2 = size(yy,1);
% 
%     Error(i) = abs(Y(size_ode,1)- yy(size2,1));
% end
% 
% polyfit(log(Error),log(h),1)
% loglog(Error,h)

%%






