[Python] Why I'm getting blank nullcline plot?

I'm trying to plot nullclines for a nonlinear system and determine fixed points from the same. However till now my code is not good enough to generate the plots. Presently I'm only able to get a blank plot only with various guess values tried for finding the root.

My objective is to generate the plots for g1 and I such that the intersecting points will be the fixed points whose coordinates are to be studied. Some coordinates as needed for the same are imported from a matlab file path included in the code.

Please advice me the corrections/fixes to be made. My code is as below

import sympy as sp
from sympy import nsolve, Abs, symbols
import numpy as np
from scipy.optimize import brentq
from pymatreader import read_mat
import matplotlib.pyplot as plt
import pandas as pd
import traceback
import math
import warnings
warnings.filterwarnings("ignore", category=RuntimeWarning)
from matplotlib.lines import Line2D
from sympy import Heaviside


# --- Load MAT file (adjust path) ---
file_path = r'E:\data&plot\data 3\set5\data1_x0_03.mat'
data = read_mat(file_path)
data1 = data['data1']
vars_ = data1['variable']
consts = data1['constant']


t = 9999
j = 99
i = 9999


# --- Pull constants ----
H0_val = consts.get('H0')
chi0_val = consts.get('chi0')
chi1_val = consts.get('chi1')
D0_val = consts.get('D0')
D1_val = consts.get('D1')
alpha_chi_val = consts.get('alphachi') # note: your const key 'alphachi'
alpha_D_val = consts.get('alphaD')
chi_val = consts.get('chi_growth', consts.get('chi_', None))
beta_val = consts.get('beta_suppression', consts.get('beta_', None))
lam_val = consts.get('lambda_suppress', consts.get('lambda', None))
mu_c_val = consts.get('critgradpressure', consts.get('critgradpressure', None))


Q_val = float(vars_['Q'][i, j])
I_val = float(vars_['turbulence'][i, j])
g1_val = float(vars_['gradpressure'][i, j])
#print(mu_c_val)
# check for missing values
missing = [k for k,v in [
('chi0', chi0_val), ('chi1', chi1_val), ('D0', D0_val), ('D1', D1_val),
('alphachi', alpha_chi_val), ('alphaD', alpha_D_val),
('chi_growth/chi_', chi_val), ('beta_suppression/beta_', beta_val),
('lambda_suppress/lambda', lam_val), ('critgradpressure', mu_c_val)
] if v is None]
if missing:
raise RuntimeError("Missing constants in MAT file: " + ", ".join(missing))

# Convert to floats (some matlab values may already be floats)
H0_val = float(H0_val)
chi0_val = float(chi0_val)
chi1_val = float(chi1_val)
D0_val = float(D0_val)
D1_val = float(D1_val)
alpha_chi_val = float(alpha_chi_val)
alpha_D_val = float(alpha_D_val)
chi_val = float(chi_val)
beta_val = float(beta_val)
lam_val = float(lam_val)
mu_c_val = float(mu_c_val)
#print(mu_c_val)
# Using v_eQ = -g1**2 => v_eQ**2 = g1**4
tot_vals = []
tot = []
I_null=[]
g1_safe=[]
def g1_dot_nullcline (g1_val, Q_val, chi0_val, chi1_val, alpha_chi_val,
I_expr, chi_val, beta_val, lam_val, mu_c_val, H0_val, tot, tot_vals):
g1_val = np.asarray(g1_val)
print(g1_val)
g1_safe = np.where(np.abs(g1_val) < 1e-12, 1e-12, g1_val)
v_eQ_sq = ( -g1_safe**2 )**2 # (-g1^2)^2 = g1^4
A = Q_val*(chi0_val + chi1_val*I_val / (1 + alpha_chi_val*(-g1_safe**2)))**(-1)
term1 = np.abs(-g1_safe - mu_c_val)
term = term1 * (term1 > 0) # Heaviside Θ
term_t = term - lam_val * v_eQ_sq
B = chi_val*term_t-beta_val*I_val**2
tot_val=A*B-H0_val
tot_vals.append(tot_val)
return (tot_vals)

for g1 in g1_safe:
# Try solving in a reasonable I-range
try:
root = brentq(nullcline_eq, -30, 30, args=(g1,))
I_null.append(root)
except:
I_null.append(np.nan)

I_null = np.array(I_null)

# --- Plot ---
plt.figure(figsize=(8,6))
plt.plot(g1_safe, I_null, lw=3)
plt.xlabel("g₁")
plt.ylabel("I (nullcline)")
plt.title("g₁-nullcline: I(g₁) from g1_dot = 0")
plt.grid(True)
plt.show()


print("total=", g1_dot_nullcline(g1_val, Q_val, chi0_val, chi1_val, alpha_chi_val, I_val,
chi_val, beta_val, lam_val, mu_c_val, H0_val, tot, tot_vals))
g1_val = np.asarray(g1_val)
print(tot_vals)

Continue reading...