% secdiff.m
% This program, secdiff.m, solves general homogeneous
%second order difference equations of the form
%     a x(n+1) + b x(n) + c x(n-1) = 0
%You are prompted to supply the constants a, b, c
% 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
a=input('input a from a x(n+1) term:    ')
b=input('input b from b x(n) term:    ')
c=input('input c from c x(n-1) term:    ')
x0=input('input initial population x(0):     ')
x1=input('input population x(1):     ')
n=input('end of time interval n:     ')
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 
      alpha=-b/(2*a);
      beta=(sqrt(-q))/(2*a);
      rho=sqrt(alpha^2+beta^2);
      theta=atan(beta/alpha);
      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



