I'm trying to write a Python code that is a numerical solver for 1-d heat conduction (using FVM) with a temperature dependent thermal conductivity. The solver has three functions I need to iterate until convergence. This is basically the algorithm from Patankar's Numerical Heat Transfer and Fluid Flow book. Conductivity function: takes an initial guess for temperature field and outputs nodal conductivities, dependent on past value of nodal temperatures.Nodal coefficients: calculates nodal coefficients based on current nodal conductivities.Nodal temperatures: uses a TDMA algorithm to calculate the nodal temperatures based on current coefficients. Once the temperatures are calculated, the values are to be cycled through to the conductivity function to continue the loop until convergence. I can make this work just fine when the conductivity is constant, but start getting issues when I try to implement this pseudocode. I get a maximum recursion depth exceeded error. My code is below:import numpy as npimport matplotlib.pyplot as pltimport copy# SOLUTION ALGORITHM# Step 1: DATA  def data(n, hw, he, TinfW, TinfE):        n = 5    hw = 10**10    he = 10**-10    TinfW = 100    TinfE = 10**12    return n, hw, he, TinfW, TinfE    # Step 2: GRID    def grid(x1, x2, data):    x1 = 0     x2 = 0.5    xe = (x2-x1) / data(n, hw, he, TinfW, TinfE)[0]    return xe# Step 3: GAMSOR     def gamsor(k, q, TDMA):    k0 = 100    beta = 20    T0 = 25       for i in range(data(n, hw, he, TinfW, TinfE)[0]):        k[i] = k0 + beta#*(TDMA(data, coeff)[i]-T0)        q[i] = 10**5    return k, q #Step 4: COEFF     def coeff(Ae, Aw, Ap, data, gamsor, grid):    for i in range(data(n, hw, he, TinfW, TinfE)[0]):               Ae[i] = Aw[i] = gamsor(k,q, TDMA(data, coeff))[0][i] /     grid(x1, x2, data)        Ap[i] = Ae[i] + Aw[i]  # Define b array - including source term and boundary conditions        b = np.zeros((data(n, hw, he, TinfW, TinfE)[0], 1))        b.fill(gamsor(k,q, TDMA(data, coeff))[1][i])            b[0] = (gamsor(k,q, TDMA(data, coeff))[1][i]*(grid(x1,x2,data)/2) +            data(n, hw, he, TinfW, TinfE)[1]*           data(n, hw, he, TinfW, TinfE)[3])        b[-1] = (gamsor(k,q, TDMA(data, coeff))[1][i]*(grid(x1,x2,data)/2)+             data(n, hw, he, TinfW, TinfE)[2]*            data(n, hw, he, TinfW, TinfE)[4])# Change values in first and last coefficients to reflect specified BC's        Aw[0] = 0        Ap[0] = Ae[0] + Aw[0] + data(n, hw, he, TinfW, TinfE)[1]        Ae[-1] = 0        Ap[-1] = Ae[0] + Aw[0] + data(n, hw, he, TinfW, TinfE)[2]    return Ae, Aw, Ap, b #Step 5: SOLVEdef TDMA(data, coeff):# Initialize T array - the solution array    T = np.zeros(data(n, hw, he, TinfW, TinfE)[0]);# Initialize recursion functions as arrays    P = np.zeros(data(n, hw, he, TinfW, TinfE)[0]);    Q = np.zeros(data(n, hw, he, TinfW, TinfE)[0]);## Step 1: evaluate at node 1    P[0] = (coeff(Ae, Aw, Ap, data, gamsor, grid)[0][0] /        coeff(Ae, Aw, Ap, data, gamsor, grid)[2][0])    Q[0] = (coeff(Ae, Aw, Ap, data, gamsor, grid)[3][0] /        coeff(Ae, Aw, Ap, data, gamsor, grid)[2][0])## Step 2: sweep from node 2 to node (n-1) (python## cell 1->(n-2) ) evaluating P and Q    for i in range(1, data(n, hw, he, TinfW, TinfE)[0]-1):        P[i] = ((coeff(Ae, Aw, Ap, data, gamsor, grid)[0][i]) /            (coeff(Ae, Aw, Ap, data, gamsor, grid)[2][i] -            coeff(Ae, Aw, Ap, data, gamsor, grid)[1][i]*P[i-1]))        Q[i] = ((coeff(Ae, Aw, Ap, data, gamsor, grid)[3][i] +            coeff(Ae, Aw, Ap, data, gamsor, grid)[1][i]*Q[i-1]) /            (coeff(Ae, Aw, Ap, data, gamsor, grid)[2][i] -            coeff(Ae, Aw, Ap, data, gamsor, grid)[1][i]*P[i-1]))## Step 3: calculate for node n    P[data(n, hw, he, TinfW, TinfE)[0]-1] = 0    Q[data(n, hw, he, TinfW, TinfE)[0]-1] = ( (coeff(Ae, Aw, Ap, data, gamsor, grid)[3][data(n, hw, he, TinfW, TinfE)[0]-1] + coeff(Ae, Aw, Ap, data, gamsor, grid)[1][data(n, hw, he, TinfW, TinfE)[0]-1]*Q[data(n, hw, he, TinfW, TinfE)[0]-2]) / (coeff(Ae, Aw, Ap, data, gamsor, grid)[2][data(n, hw, he, TinfW, TinfE)[0]-1] - coeff(Ae, Aw, Ap, data, gamsor, grid)[1][data(n, hw, he, TinfW, TinfE)[0]-1]*P[data(n, hw, he, TinfW, TinfE)[0]-2]))    T[data(n, hw, he, TinfW, TinfE)[0]-1] = (     Q[data(n, hw, he, TinfW, TinfE)[0]-1] )## Step 4: back fill, giving the temperatures    for i in range(data(n, hw, he, TinfW, TinfE)[0] - 2, -1, -1):        T[i] = P[i]*T[i+1] + Q[i]    return Tdef Temps(TDMA):    Temps = copy.copy(TDMA(data,coeff)); Temps.fill(100)    return Tempsif __name__ == '__main__': # Initialize variables    n = 0    hw = 0    he = 0    TinfW = 0    TinfE = 0    k = np.zeros(data(n, hw, he, TinfW, TinfE)[0])    q = np.zeros(data(n, hw, he, TinfW, TinfE)[0])    k0 = 0    beta = 20    T0 = 0    x1 = 0    x2 = 0    Ae = np.zeros(data(n, hw, he, TinfW, TinfE)[0])    Aw = np.zeros(data(n, hw, he, TinfW, TinfE)[0])    Ap = np.zeros(data(n, hw, he, TinfW, TinfE)[0])    b = 0    data(n, hw, he, TinfW, TinfE)    grid(x1, x2, data)    Temps(TDMA)    for _ in range(1):        gamsor(k, q)        coeff(Ae, Aw, Ap, data, gamsor, grid)        TDMA(data, coeff)#print Temps(TDMA)     print Nodal Temperature Distribution (C):, TDMA(data, coeff) 
Temperature dependent 1-d conduction in Python?
python;convergence;finite volume