How to troubleshoot numerical instability using finite difference for steady-state non-linear heat conduction equation

Question

I have a problem which I believe is numerical instability when trying to solve a heat conduction equation using finite difference. The short version is that when the parameter $I=80.3$ I get the blue curve below (positive values and with one local maximum, which is physically reasonable), but when $I=80.4$ I get the orange curve which is not physically correct. I would appreciate any help troubleshooting this problem.

Here is the long description, which I divided into physical modelling (for completeness's sake), finite difference discretization and numerical implementation (where I believe lies the problem):

1: Physical Modelling.

The problem I am trying to solve is of a thin straight metallic wire, whose both ends are held at constant temperatures $T_0$ and $T_L$. Through the wire flows a current $I$ which heats up the wire through Joule heating. The electrical resistivity is temperature-dependent as $\rho(T) = \rho_0 + \rho_1 \cdot T$

2: Finite Difference Discretization.

I divide the wire in $Nx$ nodes, for each inner node I write the following balance equation

$$q_{i-1} + q_{i+1} + q_e = 0 $$

Where $q_{i-1}$ is the heat coming from node $i-1$, $q_{i+1}$ is the heat coming from node $i+1$ and $q_e$ the heat generated by Joule heating.

I apply Fourier law $q = -k \cdot A \cdot \Delta T / \Delta x$, where $k$ is heat conductivity, $A$ is area, $\Delta T$ temperature difference and $\Delta x$ is the length between two nodes.

$$ \frac{k \cdot A}{\Delta x} \cdot (T_{i-1} - T_i) + \frac{k \cdot A}{\Delta x} \cdot (T_{i+1} - T_i) + (\rho_0 + \rho_1 \cdot T_i) \cdot \frac{\Delta x}{A} \cdot I^2 = 0 $$

Reordering I get

$$ \alpha \cdot T_{i-1} + (-2\alpha + \rho_1 \cdot \frac{\Delta x}{A} \cdot I^2) \cdot T_i + \alpha \cdot T_{i+1} = -\rho_0 \cdot \frac{\Delta x}{A} \cdot I^2 $$

3: Numerical Implementation.

I put the left side of the previous equation in a sparse tridiagonal matrix $A$, the right side in a vector $b$ and use scipy to solve the system of equations using float64. Here is the Python code:

# Imports
import numpy as np
import scipy.sparse.linalg as la
import scipy.sparse as sp
import matplotlib.pyplot as plt

#% Geometric specifications

D_L = 0.4E-3 # in m
L_L = 12E-3  # in m
T_0 = 87  # in °C
T_L = 175 # in °C

# Material constants

rho0 = 1.6E-8 # Ohm m
rho1 = 6.6E-11 # Ohm m / K
k  = 3.94E2 # W/m K 

Nx = 30 # Number of intervals

x   = np.linspace(0, L_L, Nx+1) 
d_x = x[1] - x[0]

for I in np.array( [80.3, 80.4] ):

    u = np.zeros(Nx+1) 

    # Representation of sparse matrix and right-hand side
    main  = np.zeros(Nx+1)
    lower = np.zeros(Nx)
    upper = np.zeros(Nx)
    b     = np.ones(Nx+1)

    # Construction of matrix equations

    A_c   = 3.14*(D_L)**2 /4
    alpha = k * A_c / d_x
    R     = rho0 * d_x / A_c
    Pe    = R * I**2
    b     = -Pe * b

    # Precompute sparse matrix
    main[:]  = -2*alpha + rho1 * d_x / A_c * I**2
    lower[:] = alpha
    upper[:] = alpha

    # boundary conditions
    main[0]  = 1; upper[0] = 0; b[0]  = T_0 
    main[-1] = 1; lower[-1]= 0; b[-1] = T_L 

    A = sp.diags( diagonals=[main, lower, upper], offsets=[0, -1, 1],
              shape=(Nx+1, Nx+1), format='csc' )


    u[:]   = la.spsolve(A, b)

    plt.plot( x, u, label='I='+str(I) )

plt.legend()

Does somebody know how to troubleshooting this problem? Many thanks in advance!