%
% birth death model
%
% simulates the process x-1 <- x -> x+1
%                           d(x) b(x)
% where b(x)=x
%       d(x)=x^2/xs
%
% also integrates (using very simple method) the corresponding
% deterministic process
%
% dx/dt = b(x)-d(x) (this deterministic x is denoted by xd)
%
% parameters
%
% nsteps -- number of steps
% x0 -- starting value of x
% xs -- value of x for which b(x)=d(x) -- set by input during execution
%
nsteps=1000;
x0=10;
xs=input('xs=');
xsi=1/xs;
rns=rand(2,nsteps);
rns(1,:)=log(ones(1,nsteps)-rns(1,:));
time=zeros(1,nsteps+1);x=ones(1,nsteps+1);
x(1)=x0;
timed=zeros(1,nsteps+1);xd=ones(1,nsteps+1);
xd(1)=x0;
for nt=1:nsteps
    bxd=xd(nt);dxd=xsi*xd(nt)*xd(nt);bddi=1/(bxd+dxd);
    xd(nt+1)=xd(nt)+(bxd-dxd)*bddi;timed(nt+1)=timed(nt)+bddi;
    bx=x(nt);dx=xsi*x(nt)*x(nt);bdi=1/(bx+dx);
    time(nt+1)=time(nt)-bdi*rns(1,nt);
    if x(nt)>0
        if rns(2,nt)>dx*bdi
            x(nt+1)=x(nt)+1;
        else
            x(nt+1)=x(nt)-1;
        end
    else
        x(nt+1)=0;
    end
end
plot(time,x)
hold on
plot(timed,xd,'--')
hold off


    
    