clear all
t=[0 10 20 30 40 50 60]';
N=[5.31 7.24 9.64 12.87 17.07 23.19 31.44]'
X=[ones(size(t)) t]
K=500;
c=log(N./(K-N))
d=X\c
N0=(K*exp(d(1)))/(1+exp(d(1)))
Nf=(K*N0*exp(d(2)*t))./(K-N0+N0*exp(d(2)*t))

plot(t,N,'r*',t,Nf,'-ob')
xlabel('years after 1800')
ylabel('size in millions')
legend('data','fitted model')

%prediction for 2000
%s=200;
%Ne=(K*N0*exp(d(2)*s))./(K-N0+N0*exp(d(2)*s))

%Now also add the Malthus model fitting
%p=X\log(N)
%q=exp(X*p)
%plot(t,N,'r*',t,Nf,'-ob',t,q,'-og')
%xlabel('years after 1800')
%ylabel('size in millions')
%legend('data','logistic fitted model','malthus fitted model')
%Ne2=exp(200*p(2))