% You will have to incorporate the Jacobian entries for each of the two
% fixed points (see below), i.e. modifications are required wherever
% question marks appear.
%
% S equivalent x
% I equivalent y
% xdot=x(alpha-beta*y) 
% ydot=y(beta*x-gamma*y)

global alpha beta gamma

alpha=input('enter alpha: ');
beta=0.0001;
gamma=0.03;

% first fixed point
x0=0;
y0=0;
fprintf('first fixed point (%10.5f,%10.5f) \n',x0,y0);
%--------MODIFICATIONS NEEDED HERE------------------------------------------
% include the entries of the Jacobian matrix at (x0,y0) in 
% terms of alpha, beta, gamma etc

%
% J0=[?? ?? ; ?? ??];
%---------------------------------------------------------------------------


[evector0,evalue0]=eig(J0);
fprintf('eigenvalues are %10.5f and %10.5f \n',evalue0(1,1),evalue0(2,2));
fprintf('eigenvectors are [%10.5f,%10.5f] and [%10.5f,%10.5f] \n',evector0);
% second fixed point
x1=gamma/beta;
y1=alpha/beta;
fprintf('second fixed point (%10.5f,%10.5f) \n',x1,y1);
%--------MODIFICATIONS NEEDED HERE------------------------------------------
% include the entries of the Jacobian matrix at (x0,y0) in 
% terms of alpha, beta, gamma etc

%
% J1=[?? ?? ; ?? ??];
%---------------------------------------------------------------------------

[evector1,evalue1]=eig(J1);
aa=[real(evalue1(1,1)),imag(evalue1(1,1))];
bb=[real(evalue1(2,2)),imag(evalue1(2,2))];
fprintf('eigenvalues are %9.5f + %9.5f i and %9.5f + %9.5f i \n',aa,bb);
fprintf('eigenvectors are[%10.5f,%10.5f] and [%10.5f,%10.5f] \n',evector1);


input('press enter to plot phase plane: ');

%phase plane
figure(2)
%--------------------------------------------------------------
% you may wish to MODIFY SOME OF THE FOLLOWING
%
S_0=600;
I_0=30;
R_0=0;
maxT=200;
delt=maxT/1000;
tspan=0:delt:maxT;
%
%--------------------------------------------------------------
xi=S_0;
yi=I_0;
yin=[xi;yi;R_0];
% equation for R(t) integrated also, for consistency with
% epidemic_rhs.m, but not used afterwards.
[t,yforward]=ode45('epidemic_rhs',tspan,yin);
x=yforward(:,1);
y=yforward(:,2);
plot(x,y,'g','linewidth',1);
hold on
set(gca,'fontsize',14)
xlabel('S')
ylabel('I')
title('fixed point 1: red, fixed point 2: blue')
plot(x0,y0,'ro','markersize',10)
scale=max(abs(y0),abs(y1))/5;
plot([x0-evector0(1,1)*scale x0+evector0(1,1)*scale], [y0-evector0(2,1)*scale y0+evector0(2,1)*scale],'b','linewidth',3)
  plot([x0-evector0(1,2)*scale x0+evector0(1,2)*scale], [y0-evector0(2,2)*scale y0+evector0(2,2)*scale],'b','linewidth',3)
plot(x1,y1,'bo','markersize',10)
%plot([x1-evector1(1,1)*scale x1+evector1(1,1)*scale], [y1-evector1(2,1)*scale y1+evector1(2,1)*scale],'b','linewidth',3)
%plot([x1-evector1(1,2)*scale x1+evector1(1,2)*scale], [y1-evector1(2,2)*scale y1+evector1(2,2)*scale],'b','linewidth',3)
xin=1;
ii=1;
while xin > 0 
 if(ii==1) fprintf('press enter at required coordinate location to plot trajectory: \n');end
 ii=ii+1;
 [xin,yin]=ginput(1);
 if(xin<0 | yin <0) break;end 
 [t,yforward]=ode45('epidemic_rhs',tspan,[xin;yin;R_0]);
  x=yforward(:,1);
  y=yforward(:,2);
  plot(x,y,'r')
  drawnow
  hold on
  [t,ybackward]=ode45('epidemic_rhs',-tspan,[xin;yin;R_0]);
  x=ybackward(:,1);
  y=ybackward(:,2);
  plot(x,y,'r')
end
hold off