%% Problem specifications.
% Cell format: {RHS function, interval, initial value(s), exact solution}.

% Scalar example.
ex1 = {@(t, y) -2*t*y, [0, 1], 1, @(t) exp(-t.^2)};

% System example.
ex2 = {@(t, y) [-y(2); y(1)], [0, 2*pi], [1; 0], @(t) [cos(t); sin(t)]};

ex = {ex1, ex2}; % Make a wrapper so that we can loop through the examples.

%% Run the one-step methods and collect global errors for each example.
met = {'Euler', 'RK2', 'RK4', 'exercise2'};

% Number of steps.
N = 10 * 2.^(0:10);

errors = zeros(length(ex), length(met), length(N));
for i = 1:length(ex)
    for j = 1:length(met)
        for k = 1:length(N)
            [t, y] = erk(ex{i}{1}, ex{i}{2}, ex{i}{3}, N(k), met{j});
            errors(i, j, k) = norm(ex{i}{4}(t(end)) - y(:, end), 2); % 2-norm.
        end
    end
end

%% Order of convergence-plots.
figure(1);

% Extract MATLABs colors.
Colors = get(groot, 'DefaultAxesColorOrder');

% Scalar example.
    h = diff(ex1{2}) ./ N; % Step-size.
    
    subplot(1, 2, 1)
    hold off
    loglog(h, squeeze(errors(1, :, :)), '-o')
    hold on

    % Reference lines.
    loglog(h, .05  * h,    '--', 'Color', Colors(1, :))
    loglog(h, .05  * h.^2, '--', 'Color', Colors(2, :))
    loglog(h, .005 * h.^3, '--', 'Color', Colors(4, :))
    loglog(h, .005 * h.^4, '--', 'Color', Colors(3, :))
    
    xlim([5e-5, .2])
    %ylim([1e-17, 1e-1])
    title('Scalar example', 'Interpreter', 'LaTeX', 'Fontsize', 12);
    
% System example.
    h = diff(ex2{2}) ./ N; % Step-size.
    %figure(2)
    subplot(1, 2, 2)
    hold off
    loglog(h, squeeze(errors(2, :, :)), '-o')
    hold on

    % Reference lines.
    loglog(h, 1.5   * h,    '--', 'Color', Colors(1, :))
    loglog(h,  .5   * h.^2, '--', 'Color', Colors(2, :))
    loglog(h,  .1   * h.^3, '--', 'Color', Colors(4, :))
    loglog(h,  .015 * h.^4, '--', 'Color', Colors(3, :))
    
    xlim([4e-4, 1])
    %ylim([1e-13, 15])
    title('System example', 'Interpreter', 'LaTeX', 'Fontsize', 12);

% Common plot settings.
for i = 1:length(ex)
    subplot(1, 2, i)
    xlabel(                   '$h$', 'Interpreter', 'LaTeX', 'Fontsize', 12);
    ylabel('$\| y(t_N) - y_N \|_2$', 'Interpreter', 'LaTeX', 'Fontsize', 12);
    l = legend('Euler', 'RK2', 'RK4', 'Exercise 2', ...
               '$\mathcal{O}(h)$', '$\mathcal{O}(h^2)$', '$\mathcal{O}(h^3)$', '$\mathcal{O}(h^4)$');
    set(  l, 'Color', 'none', 'Location', 'southeast', 'Interpreter', 'LaTeX', 'Fontsize', 10);
    set(gca, 'Color', 'none');
end