clear, set(0,'DefaultFigureWindowStyle','docked')
N = 40; T = 2*pi; h = T/N; A = [0 1; -1 0];
expEul = nan(2, N+1); impEul = nan(2, N+1); impMpr = nan(2, N+1);
y0 = [0; 1]; expEul(:,1) = y0; impEul(:,1) = y0; impMpr(:,1) = y0;

for ii = 1:N
    expEul(:,ii+1) = (eye(2)+h*A)*expEul(:,ii);
    impEul(:,ii+1) = (eye(2)-h*A)\impEul(:,ii);
    impMpr(:,ii+1) = (eye(2)-h*A/2)\((eye(2)+h*A/2)*impMpr(:,ii));
end
%keyboard
figure(2); subplot(1,2,1); 
plot(expEul(1,:), expEul(2,:), '*-', impEul(1,:), impEul(2,:), '*-',impMpr(1,:), impMpr(2,:), '*-', 'linewidth', 4);
axis equal
legend({'expEul', 'impEul', 'impMpr'}, 'location','northeastoutside','FontSize',24), title('y(t)', 'FontSize',24)
subplot(1,2,2), t = linspace(0, T, N+1); Q = @(y) sqrt(sum(y.^2, 1));
plot(t, Q(expEul), '*-', t, Q(impEul), '*-', t, Q(impMpr), '*-','linewidth', 4)
title('Q(y(t))','FontSize',24); fprintf('max(abs(Q(impMpr)-1)) = %e\n', max(abs(Q(impMpr)-1)))