clear
load population.dat
x= population;

N = length(x);

xn = x(1:N-1);
xnp1=x(2:N);
X=[0:max(xn)/N:max(xn)];

% plot the data
figure(1)
plot(xn,xnp1,'*',xn,xn)
xlabel('x_n','fontsize',18)
ylabel('x_{n+1}','fontsize',18)
title('raw data','fontsize',18)

% find a quadratic fit
figure(2)

[P,S]=polyfit(xn,xnp1./xn,1);
plot(xn,xnp1./xn,'r*',xn,polyval(P,xn),'b*')
xlabel('x_n','fontsize',18)
ylabel('x_{n+1}','fontsize',18)
title('quadratic fit','fontsize',18)

% find the Beverton-Holt fit
figure(3)
[P1,S]=polyfit(xn,xn./xnp1 ,1);
plot(xn,xn./xnp1,'r*',xn,polyval(P1,xn),'b*')
xlabel('x_n','fontsize',18)
ylabel('x_{n+1}','fontsize',18)
title('Beverton-Holt fit','fontsize',18)

% find the Ricker fit
figure(4)
[P2,S]=polyfit(xn,log(xnp1./xn) ,1);
plot(xn,log(xnp1./xn),'r*',xn,polyval(P2,xn),'b*')
xlabel('x_n','fontsize',18)
ylabel('x_{n+1}','fontsize',18)
title('Ricker fit','fontsize',18)

% plot all growth curves
figure(5)
plot(xn,xnp1,'*',X,polyval(P,X).*X,'r',X,X./polyval(P1,X),'b',...
    X,X.*exp(polyval(P2,X)),'g',X,X,'--')
legend('data','quadratic','B-H','Ricker','1-1')
xlabel('x_n','fontsize',18)
ylabel('x_{n+1}','fontsize',18)
title('Growth curves','fontsize',18)