function exercise5task1b

f = @(t,x) [0 -1 ; 1 0]*[x(1);-cos(x(2))];

% ODE45:
[t1,y] = ode45(f , [0,20] , [0;1] );
H1 = 0.5*y(:,1).^2 - sin(y(:,2));


% Implicit Midpoint Rule:
op = optimset('Diagnostics','off','Display','off');
h = 1/100;
t = 0;
z = [0,1];
for n=1:(20/h)
    t(n+1) = t(n)+h;
    z(n+1,:) = fsolve(@(x) x-z(n,:)-h*f(t(n)+h/2,(x+z(n,:))/2)',z(n,:),op);
end

H2 = 0.5*z(:,1).^2-sin(z(:,2));

plot(t1,H1,'b',t,H2,'r')
figure
plot(t1,y(:,1),'b',t1,y(:,2),'b--', t,z(:,1),'r',t,z(:,2),'r--')