F = @(y) [y(1)+y(2);-y(1)+y(2)];
solution = @(t) [exp(t).*cos(t);-exp(t).*sin(t)];
Y0 = [1;0];
N_v = 2.^(1:10);

M = length(N_v);

error_v = zeros(M,1);


for i = 1:M
    N = N_v(i);
    Y = Y0;
    t = 0;
    h = 1/N;
   for n = 1:N
%        Y = Y + h * F(Y);
         Y = Y + (h/2) * (F(Y) + F(Y + h * F(Y)));
   end
   error_v(i) = norm(solution(1)-Y); 
%    error_v(i) = norm(solution_approx-Y); 
    
end


loglog(1./N_v,error_v,'*-')
hold on
loglog(1./N_v,1./N_v.^(2),'--')
hold off

%%