k = .25;
A = .4;
omega = 1;
duffing = @(t, y) [y(2); A*cos(omega*t) - k*y(2) + y(1).*(1 - y(1).^2)];

% ODE45.
options = odeset('AbsTol', 1e-12, 'RelTol', 1e-12);
[t, y] = ode45(duffing, [0, 100], [0; 0], options);

figure(1); suptitle('Numerical solutions')
subplot(2, 2, 1); plot(t, y); title('ODE45')

figure(2); suptitle('Phase plots')
subplot(2, 2, 1); plot(y(:, 1), y(:, 2)); title('ODE45')

% Euler.
[t, y] = erk(duffing, [0, 100], [0; 0], 4000, 'Euler');
figure(1); subplot(2, 2, 2); plot(t, y); title('Euler')
figure(2); subplot(2, 2, 2); plot(y(1, :), y(2, :)); title('Euler')

% RK2.
[t, y] = erk(duffing, [0, 100], [0; 0], 2000, 'RK2');
figure(1); subplot(2, 2, 3); plot(t, y); title('RK2')
figure(2); subplot(2, 2, 3); plot(y(1, :), y(2, :)); title('RK2')

% RK4.
[t, y] = erk(duffing, [0, 100], [0; 0], 1000, 'RK4');
figure(1); subplot(2, 2, 4); plot(t, y); title('RK4')
figure(2); subplot(2, 2, 4); plot(y(1, :), y(2, :)); title('RK4')

% Common plot settings.
figure(1)
for i = 1:4
    subplot(2, 2, i)
    axis([-3 103 -2 2])
    xlabel('$t$', 'Interpreter', 'LaTeX', 'Fontsize', 12);
    l = legend('$u(t)$', '$u^\prime(t)$');
    set(  l, 'Color', 'none', 'Location', 'south west', 'Interpreter', 'LaTeX', 'Fontsize', 10);
    set(gca, 'Color', 'none')
end
figure(2)
for i = 1:4
   subplot(2, 2, i)
   axis([-1.6 1.6 -1.3 1.3])
   xlabel(       '$u(t)$', 'Interpreter', 'LaTeX', 'Fontsize', 12);
   ylabel('$u^\prime(t)$', 'Interpreter', 'LaTeX', 'Fontsize', 12);
   set(gca, 'Color', 'none')
end