g = 9.81;
l = 1;
A = [0,1;-g/l,0];
b = zeros(2,1);
y0 = [pi/3;0];
tstar = 8*pi*sqrt(l/g);
theta = 1/2;
mcounter = 0;
mrange = 200:200:800;
for m = mrange
  mcounter = mcounter+1;
  counter = 0;
  tstarrange = [2,4,6,8]*pi*sqrt(l/g); 
  for tstar = tstarrange
    counter = counter+1;
    y = thetametodolineare(A,b,y0,tstar,theta,m);
    err(counter,mcounter) = norm(y(:,m+1)-y0,inf);
  end
end
semilogy(mrange,10*mrange.^(-2))
t1 = text(mrange(end)+50,10*mrange(end).^(-2),'1/m^2');
hold on
i = 1;
loglog(mrange,err(i,:),'*-')
text(mrange(end)+50,err(i,end),'1 period')
for i = 2:length(tstarrange)
  loglog(mrange,err(i,:),'*-')
  text(mrange(end)+50,err(i,end),sprintf('%1d periods',i))
end
loglog(mrange,[err(1,1),err(2,2),err(3,3),err(4,4)],'r')
loglog(mrange([2,4]),[err(1,2),err(2,4)],'r')
loglog(mrange([1,2]),[err(2,1),err(4,2)],'r')
loglog(mrange,mrange/mrange(1)^2,'r')
t2 = text(mrange(end)+50,mrange(end)/mrange(1)^2,'m');
set(gca,'xtick',[200,400,600,800])
set(gca,'xticklabel',{'200','400','600','800'})
axis([150,1500,1e-6,1])
xlabel('m')
ylabel('error in infinity norm')