P0=3000; D0=100;

N=150;
P=zeros(1,N);
D=zeros(1,N);

P(1)=P0; D(1)=D0;

for n=1:N-1
    p=P(n); d=D(n);
    P(n+1)=p+0.8*(1-p/3000)*p-1.2/3000*p*d;
    D(n+1)=d-1.1*d+2.4/3000*p*d;
end;



t=1:N;
plot(t, P, 'o-g', t, D, 'x-r');