title: How to create an iterative loop of multiple Python functions?
body: <p>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. </p>
  
  <p>The solver has three functions I need to iterate until convergence: </p>
  
  <p>Conductivity function: takes an initial guess for temperature field and outputs nodal conductivities, dependent on past value of nodal temperatures.</p>
  
  <p>Nodal coefficients: calculates nodal coefficients based on current nodal conductivities.</p>
  
  <p>Nodal temperatures: uses a TDMA algorithm to calculate the nodal temperatures based on current coefficients. </p>
  
  <p>Once the temperatures are calculated, the values are to be cycled through to the conductivity function to continue the loop until convergence. </p>
  
  <p>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:</p>
  
  <pre><code>import numpy as np
  import matplotlib.pyplot as plt
  import 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: SOLVE
  def 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-&gt;(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 T
  
  def Temps(TDMA):
      Temps = copy.copy(TDMA(data,coeff)); Temps.fill(100)
      return Temps
  
  if __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) 
  </code></pre>