clear;

%% Plot eigenvalues

ks = -5:0.1:15;
lambda1 = 0.5.*(-(3-ks) + sqrt((3-ks).^2-4.*ks));
lambda2 = 0.5.*(-(3-ks) - sqrt((3-ks).^2-4.*ks));

figure(1); 
subplot(1,2,1); hold on; grid on;
plot(ks, real(lambda1));
plot(ks, real(lambda2));
xlabel('parameter k'); ylabel('imag(lambda)');
legend('lambda1', 'lambda2');

subplot(1,2,2); hold on; grid on;
plot(ks, imag(lambda1));
plot(ks, imag(lambda2));
xlabel('parameter k'); ylabel('imag(lambda)');
legend('lambda1', 'lambda2');

%% Compute the eigenvectors for k = 0.5

A = [0.5, 1; -2, -3];
[V, vals] = eig(A);

D = inv(V)*A*V;    % better: (V\A)*V

%% FitzHugh Nagumo equations

%% Plot the equilibrium as a function of I

Ve = -1.5:0.05:1.5;
Re = 1.25.*Ve + 1.5;
I = Ve.^3./3 +0.25.*Ve +1.5;

figure(3);
subplot(1,2,1);
plot(I, Ve);
xlabel('I'); ylabel('V_e');
grid on;
subplot(1,2,2); 
plot(I, Re);
grid on;
xlabel('I'); ylabel('R_e');

%% Plot the eigenvalues of the linearized system

lambda1 = 0.5*(-(-9.2+10.*Ve.^2)+sqrt((-9.2+10*Ve.^2).^2-4*(2+8*Ve.^2)));
lambda2 = 0.5*(-(-9.2+10.*Ve.^2)-sqrt((-9.2+10*Ve.^2).^2-4*(2+8*Ve.^2)));

figure(4); 
subplot(1,2,1); hold on; grid on;
plot(I, real(lambda1));
plot(I, real(lambda2));
xlabel('I'); ylabel('real(lambda)');
legend('lambda1', 'lambda2');

subplot(1,2,2); hold on; grid on;
plot(I, imag(lambda1));
plot(I, imag(lambda2));
xlabel('I'); ylabel('imag(lambda)');
legend('lambda1', 'lambda2');

%% Plot trajectories

x_dot = @(I, t, x) [10*(x(1)-(x(1)^3)/3.0-x(2)+I);
                    0.8*(-x(2)+1.25*x(1)+1.5)];

I = [0, 1.5, 2.5];

Ve = [-1.5, -0.25, 0.9];
Re = [-0.3, 1.1, 2.6];

tspan = 0:0.01:20;

figure(5);
for p=1:3
    subplot(3,1,p); hold on; grid on;
    
    x_dot_ = @(t,x) x_dot(I(p),t, x);
    
    for k=[-1, 1]
        x0 = [Ve(p)+k*0.1; Re(p)+k*0.1];
    
        [t, y] = ode45(x_dot_, tspan, x0);
        plot(t,y(:, 1), '-');
        xlabel('time t'); ylabel('voltage V');
    end
end


