p0 = 2555;
kmax = 0.026;
pmax = 12000;

ti = 0;
tf = 100;
h = 5;

% solve using ODE45
[t, pode45] = ode45(@logisticmodel,[ti:h:tf],p0,[],kmax,pmax);

% analytical solution
pop = @(t) p0*pmax./(p0 +(pmax-p0)*exp(-kmax*t)); 

% data
tdat = [1950:5:2000];
pdat = [2555 2780 3040 3346 3708 4087 4454 4850 5276 5686 6079];

% plot
figure
plot(t,pode45,t,pop(t))
hold on
plot(tdat-tdat(1),pdat,'b.')
legend('ODE45','analytical sol.','data')
xlabel('years since 1950')

% table
format short e
'year Pop(ODE45) Pop(analytical) rel.error.(%) data  residuals'
[tdat' pode45(1:length(tdat)) pop(t(1:length(tdat))) abs(pode45(1:length(tdat))-pop(t(1:length(tdat))))./pop(t(1:length(tdat))) pdat' pop(t(1:length(tdat)))-pdat']


