% integrates the FitzHugh-Nagumo model, making phase space and time evolution plots
% we use ode45, which is a Runge-Kutta method (basically a fancy Euler method)
% for Matlab reasons the ODE function itself needs to be defined at the end of the script

% DEFINE PARAMETERS
% you will have to change this
a = 0.1;

% we recommend not changing this
b=0.01;
c=0.02;

% set up start and end times of the integration
tstart = 0;
tend = 1000;

% timestep 
dt = 0.01;

% initial condition 
x_0 = [0.2 0]; 

% INTEGRATE AND PLOT
t = linspace(tstart,tend,floor((tend-tstart)/dt));

[t_out,x] = ode45(@(t,x) odefcn(t,x,a,b,c),t,x_0);

% phase space
plot(x(:,1),x(:,2));
xlabel('x')
ylabel('y')
hold on

% nullclines
x_null = linspace(-0.25,1,100);
plot(x_null,x_null.*(a-x_null).*(x_null-1));
y_null = linspace(-0.05,0.3,100);
plot(y_null.*c/b,y_null);
legend('trajectory','dx/dt=0','dy/dt=0')

figure
plot(t_out,x(:,1));
xlabel('Time')
ylabel('x')

% DEFINE ODE FUNCTION
function dxdt = odefcn(t,x,a,b,c)
    dxdt = zeros(2,1);
    dxdt(1) = -x(2)+x(1)*(a-x(1))*(x(1)-1);
    dxdt(2) = b*x(1)-c*x(2);
end

