from matplotlib.widgets import AxesWidget import numpy as np import matplotlib.pyplot as plt from scipy.integrate import solve_ivp

1. Define the parameter values

a = 0.7 b = 0.8 tau = 12.5 R = 0.1 I_ext = 0.5

1. Define the system of ODEs

def system(t, y):

   v, w = y
   dv = v - (v ** 3) / 3 - w + R * I_ext
   dw = (1 / tau) * (v + a - b * w)
   return [dv, dw]

t_span = [0, 50] initial_conditions = [(x, y) for x in np.linspace(-3, 3, 3) for y in np.linspace(-3, 3, 10)] sols = {} for ic in initial_conditions:

   sols[ic] = solve_ivp(system, t_span, ic, dense_output=True, max_step=0.01)

vmin, vmax, wmin, wmax = -2, 2, -2, 2 vs = np.linspace(vmin, vmax, 200) v_axis = np.linspace(vmin, vmax, 20) w_axis = np.linspace(wmin, wmax, 20)

v_values, w_values = np.meshgrid(v_axis, w_axis)

dv = v_values - (v_values ** 3) / 3 - w_values + R * I_ext dw = (1 / tau) * (v_values + a - b * w_values)

fig, ax = plt.subplots(figsize=(16,16))

1. integral curves

for ic in initial_conditions:

 sol = sols[ic]
 ax.plot(sol.y[0], sol.y[1], color='k', alpha=0.2, linewidth=0.5)

1. vector fields

arrow_lengths = np.sqrt(dv**2 + dw**2) alpha_values = 1 - (arrow_lengths / np.max(arrow_lengths))**0.4 ax.quiver(v_values, w_values, dv, dw, color='blue', linewidth=0.5, scale=25, alpha=alpha_values)

1. nullclines

ax.plot(vs, vs - vs**3/3 + R * I_ext, color="green", alpha=0.4, label="v nullcline") ax.plot(vs, (vs + a) / b, color="red", alpha=0.4, label="w nullcline")

ax.set_xlabel('v') ax.set_ylabel('w') ax.set_title('FitzHugh-Nagumo Model')

ax.legend() ax.set_xlim(vmin, vmax) ax.set_ylim(wmin, wmax) plt.show()