function [T]=calc_period(epsilon, h)
%
% INPUT:
%  epsilon: the small parameter
%  h:       the time step
%       Use h = 1/512 (to reduce the numerical error to below 1.0e-12)
% OUTOUT:
%  T:       the period
%
w0=[1, 0];  % Initial condition
nstep=8/h;  % Number of numerical steps
%
f_sys = @(w, t) [w(2), -sin(epsilon*w(1))/epsilon]; % Specify the ODE system to solve 
%
w=zeros(nstep+1,2);
t=zeros(nstep+1,1);
w(1,:)=w0;
%
% Use the fourth order Runge-Kutta to compute the solution
%
p=4;
d=[0,   1/2, 1/2, 1  ];
c=[0,   0,   0,   0  ;
   1/2, 0,   0,   0  ;
   0,   1/2, 0,   0  ;
   0,   0,   1,   0  ];
b=[1/6, 1/3, 1/3, 1/6];
k=zeros(p,2);
%
for j=1:nstep,
  for i=1:p,
    k(i,:)=h*f_sys(w(j,:)+c(i,1:i-1)*k(1:i-1,:), t(j)+d(i)*h);
  end
  w(j+1,:)=w(j,:)+b*k;
  t(j+1)=t(j)+h;
end
%
% Calculate the period
%
wp=w(:,2);
i0=find((wp(1:nstep)>0)&(wp(2:nstep+1)<=0));
dn=0.25/h;
T=spline(wp(i0-dn:i0+dn), t(i0-dn:i0+dn), 0);
%
%