% SySc 512: Quant Meth: Session 3, Eigenvalues and eigenvectors
% Calculation of eigenvalues for the FitzHugh-Nagumo equation
%     d/dt (x) = 10*(x - x^3/3 - y + Input)
%     d/dt (y) = (4/5)*(- y + 5*x/4 + 3/2)
%
clear
Input = 0;      %set input current (stable if Input=0, cycle if Input=2 )
% Find: isoclines
%   d/dt (x) = 0 => y = x - x^3/3 + Input
%   d/dt (y) = 0 => y = 5*x/4 + 3/2
% Graph isoclines
Xiso = -3.2:0.01:3.2;  %X for Isoclines
Isocline_x = Xiso - Xiso.^3/3 + Input;
Isocline_y = Xiso/0.8 + 1.5;

figure(3); 
Zb = plot( Xiso, Isocline_x, '-k', Xiso, Isocline_y, '--k'); 
set(Zb, 'LineWidth', 2); 
axis square;
axis([-3.5 3.5 -6 8]); xlabel('V'); ylabel('R');

% Solve isoclines simultaneously for fixed points:
%  0 = x^3/3 + 0*x^2 + x/4 + (3/2 - Input)
% Create vector of coefficients:
p = [1/3 0 1/4 3/2-Input]
% Find the roots of the polynomial:
r = roots(p)
x = r(3);
y = 5*x/4 + 3/2;
hold on                     % retain the current plot to add more than 1 function
plot(x,y,'* r')         % plot x vs t, color red
hold off                    % release hold on current plot

% Jacobian
J = [10*(1-x^2)   -10
    [      1    -4/5 ]]
% Characteristic equation:
%   (J(1,1) - lambda)*(J(2,2) - lambda) - J(1,2)*J(2,1) = 0
%    lambda^2 - (J(1,1) + J(2,2))*lambda + J(1,1)*J(2,2) - J(2,1)*J(1,2)= 0
q = [1 -(J(1,1)+J(2,2)) J(1,1)*J(2,2)-J(2,1)*J(1,2)]
% The eigenvalues are the roots oof the characteristic equation:
lambda = roots(q)
% The eigenvalues of the jacobian can be found by the Matlab command eig()
[evect eval] = eig(J)