Unable to recreate simple physics/maths paper results

[Amleto]

here is the paper:

http://www.amleto.pwp.blueyonder.co.uk/paper.pdf

I am tryng to reproduce figures 2 and 4.

my scripts they are quite simple and easy to compare with the paper. They are at the bottom of this post

I believe the author has put the inital condition sin and cos the wrong way around, but even so, my scripts wont work for delta <0.2 really (he uses as low as 0.01)

can anyone see an error? or a better method?

Thanks

=======

http://www.amleto.pwp.blueyonder.co.uk/help1.m

% Tom Collings 2004

global beta alpha gamma


% Solver setup----------------

tol=1e-6;

t0 = 0;
n = 12; % the number of cycles

tf = 2*pi*n;
stepsize = 2*pi/250; %please keep to fractions of 2*pi
tspan = 0:stepsize:tf;


delta=0.3;
theta0 = 20 * pi/180;
rho0 = (1+delta)*sin(theta0);
z0 = (1+delta)*cos(theta0);


% R == a in the paper
initcons = [rho0 z0];


options = odeset('reltol',tol);

alpha=0.5;
beta=0.2;
gamma = pi/3;


[t,r] = ode45('help1func',tspan,initcons,options);

% ==================================================
% plotting  =========================================
% ====================================================



% Particle Path ===================
% ===============================

rmax = 1+alpha;
rmin = 1-alpha;

figure
plot(r(:,1),r(:,2))

hold;

tau=linspace(0,2*pi,50);
plot(rmin*cos(tau),rmin*sin(tau) - 0,'k-.')
plot(rmax*cos(tau),rmax*sin(tau) - 0,'k-.')

axis equal;


beep

ii=find( t/2/pi==fix(t/2/pi));

tt = t(ii);

figure
hold on
plot(r(ii,1),r(ii,2),'o');
plot(rmin*cos(tau),rmin*sin(tau) - 0,'k-.');
plot(rmax*cos(tau),rmax*sin(tau) - 0,'k-.');
axis equal

http://www.amleto.pwp.blueyonder.co.uk/help1func.m

function dr = help1func(t,r)

global beta alpha gamma


h = beta * sin(t + gamma);
hdot = beta* cos(t + gamma);


%rc == r in the paper
rc = sqrt( r(1).^2 + (r(2) - h).^2 );


rad=1+alpha*sin(t);
raddot= alpha*cos(t);


% r(1) == rho
% r(2) == z
dr(1) =  raddot.*rad.^2.*r(1)./rc.^3        + 3/2* hdot.*rad.^3.* r(1).*( r(2) - h )./rc.^5;
dr(2) =  raddot.*rad.^2*( r(2) - h )./rc.^3 + 1/2* hdot.*rad.^3.* ( 3*( r(2) - h ).^2./rc.^5 - 1./rc.^3 );


     
dr=dr';

%tleft = (tf-t)/tf*100