FITZHUGH_NAGUMO

=FITZHUGH_NAGUMO(v_initial, w_initial, a, b, tau, i_ext, t_start, t_end, [timesteps], [solve_ivp_method])

=FITZHUGH_NAGUMO(0.0, 0.0, 0.7, 0.8, 12.5, 0.5, 0.0, 1.0, 10)

=FITZHUGH_NAGUMO(0.0, 0.0, 0.7, 0.8, 12.5, 0.5, 0.0, 1.0, 10, "RK23")

=FITZHUGH_NAGUMO(1.0, 0.5, 0.7, 0.8, 12.5, 0.5, 0.0, 2.0, 10, "RK45")

=FITZHUGH_NAGUMO(-1.0, 0.2, 0.5, 0.7, 10.0, 0.3, 0.0, 1.5, 10)

from scipy import integrate as scipy_integrate
import numpy as np
 
 
def fitzhugh_nagumo(v_initial, w_initial, a, b, tau, i_ext, t_start, t_end, timesteps=10, solve_ivp_method='RK45'):
    """
    Wrapper around scipy.integrate.solve_ivp that numerically solves the FitzHugh-Nagumo system of
    ordinary differential equations for neuron action potentials.
 
    This function wraps the scipy.integrate.solve_ivp solver. See:
    https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html
 
    Args:
        v_initial: Initial membrane potential.
        w_initial: Initial recovery variable.
        a: Parameter a.
        b: Parameter b.
        tau: Time scale parameter.
        i_ext: External current.
        t_start: Start time of integration.
        t_end: End time of integration.
        timesteps: Number of timesteps to solve for. Default is 10.
        solve_ivp_method: Integration method ('RK45', 'RK23', 'DOP853', 'Radau', 'BDF', 'LSODA'). Default is 'RK45'.
 
    Returns:
        2D list with header row: t, V, W. Each row contains time and variable values, or an error message (str) if input is invalid.
 
    This example function is provided as-is without any representation of accuracy.
    """
    # Validate input types
    try:
        v0 = float(v_initial)
        w0 = float(w_initial)
        a = float(a)
        b = float(b)
        tau = float(tau)
        i_ext = float(i_ext)
        t0 = float(t_start)
        t1 = float(t_end)
        ntp = int(timesteps)
    except Exception:
        return "Invalid input: All initial values, parameters, and timesteps must be numbers."
    if t1 <= t0:
        return "Invalid input: t_end must be greater than t_start."
    if tau <= 0:
        return "Invalid input: tau must be positive."
    if ntp <= 0:
        return "Invalid input: timesteps must be a positive integer."
    if solve_ivp_method not in ['RK45', 'RK23', 'DOP853', 'Radau', 'BDF', 'LSODA']:
        return "Invalid input: solve_ivp_method must be one of 'RK45', 'RK23', 'DOP853', 'Radau', 'BDF', 'LSODA'."
    # Create time array for evaluation
    t_eval = np.linspace(t0, t1, ntp)
    
    # FitzHugh-Nagumo ODE system
    def fhn_ode(t, y):
        v, w = y
        dv_dt = v - v**3 / 3 - w + i_ext
        dw_dt = (v + a - b * w) / tau
        return [dv_dt, dw_dt]
    # Integrate
    try:
        sol = scipy_integrate.solve_ivp(
            fhn_ode,
            [t0, t1],
            [v0, w0],
            method=solve_ivp_method,
            dense_output=False,
            t_eval=t_eval
        )
    except Exception as e:
        return f"Integration error: {e}"
    if not sol.success:
        return f"Integration failed: {sol.message}"
    # Format output: header row then data rows
    result = [["t", "V", "W"]]
    for i in range(len(sol.t)):
        t_val = float(sol.t[i])
        v_val = float(sol.y[0][i])
        w_val = float(sol.y[1][i])
        # Disallow nan/inf
        if any([
            not isinstance(t_val, float) or not isinstance(v_val, float) or not isinstance(w_val, float),
            t_val != t_val or v_val != v_val or w_val != w_val,  # NaN check
            abs(t_val) == float('inf') or abs(v_val) == float('inf') or abs(w_val) == float('inf')
        ]):
            return "Invalid output: nan or inf detected."
        result.append([t_val, v_val, w_val])
    return result