%
%                     Fitzhugh-Nagumo Model
%                       Hopf Bifurcations
%                Victor Matveev, Oct 16 2015
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure('Position',[1 200 700 700]); % Specify window size
A = 1/2; C = 1/2;
Icr = -sqrt(1 - A*C)*(1/C - 2/3 - A*C/3);
tTotal = 300;
L = 2.5; LW = 3;
fs = 18;

for I = [-2:0.02:2];  % Loop over different values of applied current I
    
        f = @(t,y) [ y(1) - y(1).^3/3 - y(2) + I; A*(y(1) - C*y(2)) ];
        g = @(y) f(0,y);

        fp = fsolve(g,[0 0]);       % Find the fixed point
        Vss = fp(1); Wss = fp(2);   % Get the steady-state V and W values from "fp"
        J = [ [1 - Vss^2, -1]; [A, -A*C]];  % The Jacobian
        [Vectors, Lambda] = eig(J);                          % Eigenvalues of Jacobian

        V = linspace(-L, L, 60);
        WN = V / C;  % W-nullcline
        VN = V - V.^3/3 + I;   % V-nullcline

        %%%%%%%%%%%%%%%%%%%%%%%     PANEL  1     %%%%%%%%%%%%%%%%%%%%%%%%

        subplot(2,2,1); 
        hold off;
        plot(V, WN, '-', 'linewidth', 2, 'color', [0 0.5 0]);    % plot the W-nullcline
        hold on;
        plot(V, VN, '-', 'linewidth', 2, 'color', [0 0 0.2]);   % plot the V-nullcline

        IC = [-1 0];
        IC = [Vss+0.2, Wss];
        IC = [-0.6, -2.57];
        [T Y] = ode45(f, [0 tTotal], IC);
        plot(Y(:,1), Y(:,2), 'm-', 'linewidth', 1.5);

        set(gca, 'fontsize', fs);
        title(['Phase plane, I=', num2str(I)], 'fontsize', fs);
        xlabel('V(t)'); ylabel('W(t)');
        axis([-L L -LW LW]);

        vv = linspace(-L, L, 8);
        ww = linspace(-LW, LW, 8);
        [VV, WW] = meshgrid(vv, ww);
        quiver(VV, WW, VV-VV.^3/3-WW+I, A*(VV-C*WW));
        
        %%%%%%%%%%%%%%%%%%%%%%%     PANEL  2     %%%%%%%%%%%%%%%%%%%%%%%%

        subplot(2,2,2); hold off; 
        plot(T, Y(:,1), 'm-', 'linewidth', 1.5); hold on;
        plot([0 tTotal],[Vss Vss],'k:', 'linewidth', 1.5);
        set(gca, 'fontsize', fs);
        title('Voltage vs time');
        xlabel('time'); ylabel('V(t)');
        axis tight;

        %%%%%%%%%%%%%%%%%%%%%%%     PANEL  3     %%%%%%%%%%%%%%%%%%%%%%%%

        subplot(2,2,3); hold on; 
        if real(Lambda(1,1)) > 0
            color = [1 0 0];   % Red color  ([R G B] values)
        else 
            color = [0 0 1];   % Blue color ([R G B] values)
        end;
        plot(real(Lambda(1,1)),imag(Lambda(1,1)), 'o', 'color', color, 'markersize', 4);
        plot(real(Lambda(2,2)),imag(Lambda(2,2)), 'o', 'color', color, 'markersize', 4);
        set(gca, 'fontsize', fs);
        title(['Real(\lambda)=',num2str(real(Lambda(1,1)))], 'fontsize', fs);
        xlabel('Re(\lambda)'); ylabel('Im(\lambda)');
        axis tight;

        %%%%%%%%%%%%%%%%%%%%%%%     PANEL  4     %%%%%%%%%%%%%%%%%%%%%%%%

        subplot(2,2,4); hold on; set(gca, 'fontsize', fs);
        [T Y] = ode45(f, [0 tTotal], Y(length(T),:) );  % Run some more in order
        [T Y] = ode45(f, [0 tTotal], Y(length(T),:) );  % to settle at equilibrium
        plot(I, min(Y(:,1)), 'mo', 'markersize', 4);       % minimum of V(t): a magenta point
        plot(I, max(Y(:,1)), 'mo', 'markersize', 4);       % maximum of V(t): a magenta point
        plot(I,Vss,'o','color',color, 'markersize', 4);    % Fixed point: a red or blue point
        title('Bifurcation diagram (V_{min}, V_{max}, V*)', 'fontsize', 12);
        xlabel('I'); ylabel('V equilibrium');
        axis tight;
        drawnow;

end;



