% This is the main function to solve 
% the Fitzhugh-Nagumo model 


global a epsi gamma

a = 0.3;
epsi = 0.01;
gamma = 2.5;



options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4]);

t0 = 0;
tfinal = 100;


v_ini = [0.4 0.5 0.6];
w_ini = [0.0 0.0 0.0];

for ind = 1:3

y_ini = [v_ini(ind) w_ini(ind)];

[T,y] = ode45(@FitzhughNagumo_eqn,[t0 tfinal],y_ini,options);

% plot N2 versus N1

figure(1);hold on; plot(y(:,1),y(:,2));
xlabel('x_1')
ylabel('x_2')
axis equal; box on

figure(2); hold on;
plot(T,y(:,1),T,y(:,2))
legend('v','w')

end


% plot the nullcline

v = -0.1:0.01:1.05
w1 = -v.*(v-a).*(v-1);
figure(1);hold on; plot(v,w1,'k')
w2 = v/gamma;
figure(1);hold on; plot(v,w2,'k')