
[trk,xrk] = rk_classical(@linsystem, 0, 0.1, [2 1], .0025);
[t1,x1] = ode12_1(@linsystem, 0, 0.1, [2 1]);
[t2,x2] = ode12_2(@linsystem, 0, 0.1, [2 1]);
[t,x] = ode12(@linsystem, 0, 0.1, [2 1]);

plot(trk, xrk(:,1), '-kx', 'DisplayName', 'Runge-Kutta');
hold on;
plot(t1, x1(:,1), '-mx', 'DisplayName', 'ode12\_1');
plot(t2, x2(:,1), '-bo', 'DisplayName', 'ode12\_2');
plot(t, x(:,1), '-rs', 'DisplayName', 'ode12');
legend('Location', 'east');
xlabel('t');
ylabel('x_1');
ylim([1.5 6.5]);

figure;
plot(trk, xrk(:,2), '-kx', 'DisplayName', 'Runge-Kutta');
hold on;
plot(t1, x1(:,2), '-mx', 'DisplayName', 'ode12\_1');
plot(t2, x2(:,2), '-bo', 'DisplayName', 'ode12\_2');
plot(t, x(:,2), '-rs', 'DisplayName', 'ode12');
legend('Location', 'east');
xlabel('t');
ylabel('x_2');
ylim([-4 2]);

disp('Number of time points:')
fprintf('  RK      = %d\n', length(trk));
fprintf('  ode12_1 = %d\n', length(t1));
fprintf('  ode12_2 = %d\n', length(t2));
fprintf('  ode12   = %d\n', length(t));
