ode systems including linear and nonlinear equations

Question

I am dealing with ode equations that are both nonlinear and linear and change over time. I have solved them using scipy integrate odint, but with changing t_span the result will change unexpectedly e.g. results at time span of between zero to one will not correspond to time span zero to 0.

import numpy as np
import math
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# Define parameters and functions
t0 = 0.0
T = 1
N = 10 # Number of time steps
A = np.array([10**14, 4 * 10**13, 10**15, 10**12])
E = np.array([2.24e-19, 2.03e-19, 2.24e-19, 1.4e-19])
h = np.array([1.714e6, 257e5, 3.14e5, 3.78e4])
kb = 1.38 * (10**(-23))
z_t0 = 0.33
C = 1.1  # Battery heat capacity (J/K)
R_conv = 10  # Convective heat transfer coefficient (W/K)
R_rad = 0.8  # Radiative heat transfer coefficient (W/K)
M = 0.72

def ode_system(y, t):
    x_a, x_c, x_s, z, soc, tem = y
    # if x_a <= 0 or x_c <= 0 or x_s <= 0 or soc <= 0:
    #     return [0, 0, 0, 0, 0, 0]
    dx_adt = -A[0] * x_a * math.exp(-E[0] / (tem * kb))
    dx_cdt = A[1] * x_c * (1 - x_c) * math.exp(-E[1] / (tem * kb))

    try:
        deli1 = math.exp(-E[2] * (-z) / (tem * kb * z_t0))
    except OverflowError:
        deli1 = np.inf

    dx_sdt = -A[2] * x_s * deli1
    deli2 = np.exp(-E[2] * (-z) / (tem * kb * z_t0))
    try:
        deli2 = math.exp(-E[2] * (-z) / (tem * kb * z_t0))
    except OverflowError:
        deli2 = np.inf
    dzdt = A[0] * x_s * deli2
    dsocdt = -A[3] * (1 - x_c) * x_a * math.exp(-E[3] / (tem * kb))

    def total_heat_generation():
        Q_xa = 0.13 * h[0] * dx_adt
        Q_s = -0.13 * h[1] * dx_sdt
        Q_c = 0.29 * h[2] * dx_cdt
        Q_isc = -0.37 * h[3] * dsocdt
        Q_TR = Q_xa + Q_s + Q_c + Q_isc
        return Q_TR

    dTdt = total_heat_generation() / M * C + (R_conv * (tem - 298) - R_rad * (tem ** 4 - 298 ** 4)) / M * C*100/32

    return [dx_adt, dx_cdt, dx_sdt, dzdt, dsocdt, dTdt]

# Set the time span, initial conditions
t_span = np.linspace(t0, T, N+1)
initial_conditions = [0.75, 0.04, 0.015, 0.033, 1, 298]

# Solve the ODEs using odeint
y_values = odeint(ode_system, initial_conditions, t_span)

# Extract individual solutions
x_av, x_cv, x_sv, z_v, socv, temv = y_values.T

# Example: Plotting the results
plt.figure(figsize=(10, 6))
plt.subplot(2, 1, 1)
plt.plot(t_span, x_av, label='x_a')
plt.plot(t_span, x_cv, label='x_c')
plt.plot(t_span, x_sv, label='x_s')
plt.plot(t_span, socv, label='SOC')
plt.legend()
plt.title('Concentration Profiles')
plt.xlabel('Time')
plt.ylabel('Concentration')
plt.show()
print("x_av[1]:", x_av[1])
print("socv[1]:", socv[1])
x_av