%
%  code31.m -- gumdrop period (28 nov 2001 - djm)
%
%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

magic = 0.09814515000782;

%  ODE file, IVs, final time, method & options

Xode = 'code29a';
method = 'ode45';  tol = 1e-12;

Tbl = [];  Tend = [];  Tmg = [];

Ups = 2:2:50;

for j=1:size(Ups,2)
  disp([num2str(Ups(j))]);
  Tbl = [Tbl ; 2*sqrt(Ups(j)) + (1+log(4))/2/sqrt(Ups(j))];  
  Tmg = [Tmg ; Tbl(j) + 2*magic/sqrt(Ups(j)^3)];  
  IVs = [0 Ups(j)];

  options = odeset('reltol',tol,'abstol',tol,'refine',1, ...
		'stats','off','events','on');

%  ODE solve with period & amplitude determination

  [t y] = feval(method,Xode,[0 2*Tbl(end)],IVs,options);
  Tend = [Tend ; t(end)];
end

figure(4);  clf
plot(log(sqrt(Ups)),log(2*abs(Tbl-Tend)),'x')
hold on
plot(log(sqrt(Ups)),log(2*abs(Tmg-Tend)),'o')
plot([0.8 2.8],[-3  -9],'r')
plot([0.8 2.8],[-5 -15],'r')

title('\bf error scaling (P = period)')
xlabel('\bf log(A)')
ylabel('\bf log|P_{rk}-P_{bl}|')
text(1.3,-4,'\bf slope = -3')
text(1.5,-8,'\bf slope = -5')

axis([0.2 2 -12 -2])