% Creates a table plotting Effort (function evaluations) vs. Error

% Function and initialization info
dfun = 'orbit'; 
tspan = [0 6.19216933];     % must be exactly 1 period
y0 = [1.2; 0; 0; -1.04935750983031990726];

% Error parameters
ErrorHigh = 9.5;
ErrorLow = 5;
ErrorStep = 0.2;            % must be non-zero

% RK(4,4) Formulas
% op = Our formula (c2 = 0.3578, c3 = 0.5915)
% kt = Kutta formula (c2 = 1/3, c3 = 2/3)
% kn = Kuntzmann formula (c2 = 2/5, c3 = 3/5)
% ra = Ralston formula (c2 = 0.4, c3 = 0.45573725)
% cl = Classic formula (c2 = 1/2, c3 = 1/2)
% hj = Hull/Johnston formula (c2 = 0.35, c3 = 0.45)
% od = ode45

for k = ErrorLow:ErrorStep:ErrorHigh
    a = (k - ErrorLow)/ ErrorStep + 1;  % a is an integer starting at 1
    a = int8(a);                       % that increases at increments of 1
    
    % RK44Auto_3 is an edited version of RK44Auto that 
    % outputs the number of deriv evaluations and the error
    % Stats output has also been disabled
                                        
    [ef_op(a),er_op(a)] = RK44Auto_3(dfun, tspan, y0, 0.3578, 0.5915, 10^(-k), 10^(-k));
    
    [ef_kt(a),er_kt(a)] = RK44Auto_3(dfun, tspan, y0, 1/3, 2/3, 10^(-k), 10^(-k));
    
    [ef_kn(a),er_kn(a)] = RK44Auto_3(dfun, tspan, y0, 2/5, 3/5, 10^(-k), 10^(-k));
    
    [ef_ra(a),er_ra(a)] = RK44Auto_3(dfun, tspan, y0, 0.4, 0.45573725, 10^(-k), 10^(-k));
    
    [ef_cl(a),er_cl(a)] = RK44Auto_3(dfun, tspan, y0, 1/2, 1/2, 10^(-k), 10^(-k));
    
    [ef_hj(a),er_hj(a)] = RK44Auto_3(dfun, tspan, y0, 0.35, 0.45, 10^(-k), 10^(-k));

    
    % ode45
    options = odeset('AbsTol',10^(-k),'RelTol',10^(-k));
    [x_od,y_od,stats] = ode45(dfun, tspan, y0, options);
    er_od(a) = -log10(norm([y_od(length(y_od),:) - y_od(1,:)]));
    ef_od(a) = stats(3);
    
end

plot(er_op,ef_op,'b',er_kt,ef_kt,'r',er_kn,ef_kn,'g',er_ra,ef_ra,'c',...
    er_cl,ef_cl,'m',er_hj,ef_hj,'y',er_od,ef_od,'k')
xlabel('Error')
ylabel('Effort')
h = legend('REU06','Kutta','Kuntzmann','Ralston','Classic','Hull/Johnston','ode45','Location','NorthWest');
set(h);