clear all
theta0 = pi/3;
g = 1;
l = 1;
tstar = 3*pi*sqrt(l/g);;
timesteps = 100;
h = tstar/timesteps;
t = 0;
y(:,1) = [theta0,0];
A = [0,1;-g/l,0];
theta = 1/2;
for n = 1:timesteps
  y(:,n+1) = (eye(2)-(1-theta)*h*A)\((eye(2)+theta*h*A)*y(:,n));
%  plot(theta(1,n+1),theta(2,n+1),'*')
%  axis([-2,2,-2,2])
%  pause(0.05)
end
plot(y(1,:),y(2,:))
axis([-2,2,-2,2])
axis square
xlabel('y1')
ylabel('y2')
