function sol1

figure(1)
clf
ts = ([0:0.01:0.99 0.999])*pi;
T = zeros(size(ts));
TE = zeros(size(ts));
T(1) = 2*pi;
for j = 2:length(ts)
  T(j) = 4*integral(@(t) 1./sqrt(2*(cos(t)-cos(ts(j)))),0,ts(j));
end
% matlab uses the parameter
[K,E] = ellipke(sin(0.5*ts).^2);
Te = 4*K;
figure(1)
plot(ts/pi,T/(2*pi),ts/pi,Te/(2*pi))
axis([0 1 0 4])
xlabel('\theta_*/\pi')
ylabel('(l/g)^{1/2} T/(2\pi)')


figure(3)
clf
hold on
rr = 2.4:0.01:4;
plot(rr,1-1./rr,'r','linewidth',2)
r2 = 3:.01:4;
sp = 0.5 + 0.5./r2 + 0.5*sqrt(1-2./r2-3./r2.^2);
sm = 0.5 + 0.5./r2 - 0.5*sqrt(1-2./r2-3./r2.^2);
plot(r2,sp,'b','linewidth',2)
plot(r2,sm,'b','linewidth',2)
x = zeros(200,1);
for r = rr
  x(1) = 0.5;
  for j = 1:500
    x(j+1) = r*x(j)*(1-x(j));
  end
  plot(r,x(101:200),'.k')
end
hold off
xlabel('r')
ylabel('x')


figure(4)
clf
x0 = -2:0.25:3;
y0 = -1:0.25:3;
hold on
options = odeset('Events',@f4e);
warning off
for j = 1:length(x0)
  for k = 1:length(y0)
    [t,y]=ode45(@(t,z) f4(t,z),[0 10],[x0(j) y0(k)],options);
    plot(y(:,1),y(:,2),'b')
    [t,y]=ode45(@(t,z) f4(t,z),[0 -10],[x0(j) y0(k)],options);
    plot(y(:,1),y(:,2),'b')
  end    
end
warning on
plot([0 2 0 -1],[0 2 1 2],'rx')
axis equal
axis([-1.5 2.5 -0.5 2.5])
xlabel('x')
ylabel('y')
hold off
% alternative approach
%N = chebop(0,20);
%N.op = @(t,x,y) [diff(x)-(x-y)*(1-x-y)/3 ; diff(y)-x*(2-y)]
%quiver(N,[-0.5 3.5 -0.5 3.5],'normalize',1)


function zdot = f4(t,z)

x = z(1); y = z(2);
zdot = [(x-y)*(1-x-y)/3 ; x*(2-y)];


function [value,isterminal,direction] = f4e(t,z)

x = z(1); y = z(2);
value = [x+1.5 x-2.5 y+0.5 y-2.5];
isterminal = [1 1 1 1];
direction = [-1 1 -1 1];

