function dydt=fhn(t,y,opt)
% dydt=fhn(t,y,opt)
% y=[V;R];
% opt.u = input
% opt.a, opt.b and opt.c are parameters
%
% FitzHugh-Nagumo neural spike potential equations
% dV/dt = c(V-V^3/3+R)+u(t)
% dR/dt = -1/c(V-a+b*R)
%
% S.Jbabdi
a = opt.a;
b = opt.b;
c = opt.c;
u = opt.u;

% nearest neighbour interpolation
[~,i] = min(abs(opt.tint-t));
ut = u(i);

dydt = 0*y;
dydt(1) = c*(y(1)-y(1)^3/3+y(2)) + ut;
dydt(2) = -1/c*(y(1)-a+b*y(2));
