The warnings you’re seeing (RuntimeWarning: overflow encountered in scalar add/divide, invalid value encountered in scalar add/subtract) and the plot blowing up to ±1e6 or worse (with wild oscillations and negative values) are classic signs of numerical instability in the explicit finite difference scheme.

# Force diffusion stability instead of pure advection CFL
diff_number_max = 0.45
dt_diff = diff_number_max * dx**2 / D
dt = min(CFL * dx / U, dt_diff)          # take the more restrictive one

import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse import diags
from scipy.sparse.linalg import spsolve
# ─── Parameters ───────────────────────────────────────────────
Lx = 100.0
Nx = 601               # bit coarser to make matrix smaller/faster
dx = Lx / (Nx - 1)
x = np.linspace(0, Lx, Nx)
U   = 2.0
D   = 0.5
λ   = 0.12
α   = 3.0
T_end = 40.0
CFL = 0.8
dt  = CFL * dx / U      # advection controls dt — diffusion is now unconditional
Nt  = int(T_end / dt) + 1
print(f"dx = {dx:.4f},  dt = {dt:.4f},  Nt = {Nt}")
def S(t):
    return np.exp(-((t - 8.0)/2.5)**2)
# Initial condition
C = np.zeros(Nx)
C[x <= 10] = 1.0
# ─── Build tridiagonal matrix for Crank-Nicolson diffusion ────
r = D * dt / (2 * dx**2)           # half because CN averages two levels
main_diag  = np.ones(Nx) * (1 + 2*r + λ*dt)
off_diag   = np.ones(Nx-1) * (-r)
# We'll modify boundaries later
A = diags([off_diag, main_diag, off_diag], [-1, 0, 1], format='csr')
C_history = [C.copy()]
times = [0.0]
# ─── Time loop ─────────────────────────────────────────────────
for n in range(Nt):
    t = (n+1) * dt
    # Right-hand side: explicit part (advection + old diffusion + decay + source)
    b = np.zeros(Nx)
    # Upwind advection (explicit)
    adv = np.zeros(Nx)
    adv[1:] -= U * (C[1:] - C[:-1]) / dx
    # Explicit half diffusion + decay + source
    diff_old = D * (np.roll(C, -1) - 2*C + np.roll(C, 1)) / dx**2
    diff_old[0] = diff_old[-1] = 0  # boundaries
    b = C + dt * (adv + 0.5 * diff_old - λ * C + α * S(t))
    # Implicit half diffusion is already in A
    # Boundaries: Dirichlet left = 0, Neumann right ≈ zero gradient
    A = A.copy()  # reset if needed
    b[0]  = 0.0
    A[0, :] = 0
    A[0, 0] = 1
    A[-1, :] = 0
    A[-1, -2] = -1
    A[-1, -1] = 1
    b[-1] = 0   # dC/dx ≈ 0 → C[-1] = C[-2]
    # Solve
    C_new = spsolve(A, b)
    C = C_new
    if n % max(1, int(1.0/dt)) == 0:   # save ~every 1 time unit
        C_history.append(C.copy())
        times.append(t)
# ─── Plot ──────────────────────────────────────────────────────
fig, ax = plt.subplots(figsize=(11, 5))
for C_snap, ts in zip(C_history[::max(1,len(C_history)//12)], times[::max(1,len(times)//12)]):
    ax.plot(x, C_snap, lw=1.2, label=f"t={ts:.1f}")
ax.set_xlabel("x")
ax.set_ylabel("C(x,t)")
ax.set_title(f"Upwind adv + CN diff  |  U={U}  D={D}  λ={λ}  α={α}")
ax.grid(True, alpha=0.25)
ax.legend(ncol=2, fontsize=9)
plt.show()