Simulation of the FitzHugh-Nagumo model

The FitzHugh-Nagumo model is an example of a two-dimensional dynamical system that is simple and behaves in essence like the Hodgkin-Huxley model regarding to fast-slow phase plane. So it has two variables: one fast excitation variable $v$ and one slow recovery variable $w$. The system of equations is:

$$\begin{aligned}\frac{{\dd v}}{{\dd t}} &= f(v) - w + {I_{{\rm{stim}}}}\\[0.25cm] \frac{{\dd w}}{{\dd t}} &= \varepsilon \cdot (v + \beta - \gamma \cdot w)\end{aligned}$$ where $\varepsilon$ is a small parameter and $f$ is a suitable function.

Traditional choices are:

$$f(v) = v - {\tfrac{1}{3}}{v^3},\quad \varepsilon = 0.08,\quad \beta = 0.7,\quad \gamma = 0.8$$