clear all
global Pi alpha mu kappa beta gamma muW

Pi=37;
alpha=1;
mu=0.0142;
kappa=365*24;
beta=0.0255;
gamma=100000;
muW=26;
r=0.9;

t0=0;
tq=[];
yq=[];
tfinal=12;
x0=[Pi/mu 0 0 200 200];
tq=[];
xq=[];
mq=[];
t0=0;
reps = 6;
tau=1;

for i=1:reps
 tspan=[t0 t0+tau];
 [t,x] = ode45(@GWDf,tspan,x0);
 n=length(x);
 x0=x(n,:);
 x0(4)=(1-r).*x(n,4);
 x0(5)=(1-r).*x(n,5);
 tf=t(n);
 tq=[tq;t(1:n)];
 xq=[xq;x((1:n),:)];
 mq=[mq;mean(x)];
 t0=tf;
end
%plot(tq,xq(:,5))
plot(tq,xq(:,1),tq,xq(:,2),tq,xq(:,3),tq,xq(:,4))

xlabel('time (years)')
ylabel('Population')
