% seedmod.m
% This program, seedmod.m, solves the seedbank model
% second order difference equation from chapter 1 of
% Edelstein-Keshet (Math Models in Biology)
% The equation is of the general form
%     a x(n+1) + b x(n) + c x(n-1) = 0
% where for this case a=1 and b=-alpha*sigma*gamma
% and c=-beta*(sigma^2)*gamma*(1-alpha)
% You are prompted to supply the constants 
% alpha, beta, gamma, sigma
% as well as x(0), x(1) and a time interval over
% which to plot the solution. This program handles
% all cases for the roots of the characteristic
% equation.
!c:
s=1;
while s>0
alpha=input('input alpha=1 year old germination fraction:    ')
beta=input('input beta=2 year old germination fraction:    ')
gamma=input('input gamma=number seeds per plant:    ')
sigma=input('input sigma=overwinter seed survival fraction:    ')
x0=input('input initial population x(0):     ')
x1=input('input population x(1):     ')
n=input('end of time interval n:     ')
a=1.;
b=-alpha*sigma*gamma
c=-beta*(sigma^2.)*gamma*(1.-alpha)
x=zeros(n+1,1);
t=zeros(n+1,1);
x(1)=x0;
x(2)=x1;
t(2)=1;
q=b*b-4.*a*c;
if q>0 
    r1=(-b+sqrt(q))/(2*a);
    r2=(-b-sqrt(q))/(2*a);
    c2=(x1+x0*r1)/(r1+r2)
    c1=c2-x0
      for i=2:n
      t(i)=i;
      x(i+1)=c1*r1^i+c2*r2^i;
      end
     t(n+1)=n;
     plot(t,x,t,x,'o'),pause
end
if q==0 
     r1=-b/(2*a);
     c1=x0;
     c2=-x0+x1/r1;
      for i=2:n
      t(i)=i;
      x(i+1)=(c1+c2*i)*r1^i;
      end
     t(n+1)=n;
     plot(t,x,t,x,'o'),pause
end
if q<0 
      al=-b/(2*a);
      be=(sqrt(-q))/(2*a);
      rho=sqrt(al^2+be^2);
      theta=atan(be/al);
      c1=x0;
      c2=((x1/rho)-x0*cos(theta))/sin(theta);
      for i=2:n
      t(i)=i;
      x(i+1)=(c1*cos(i*theta)+c2*sin(i*theta))*rho^i;
      end
     t(n+1)=n;
     plot(t,x,t,x,'o'),pause
end
s=input('To stop type a 0 - else type a 1'),pause
end