title: Solving non-linear coupled differential equations in python
body: <p>I am working on simulation of a system that contains coupled differential equations. My main aim is to solve the mass balance in steady condition and feed the solution of steady state as initial guess for the dynamic simulation.
  There are basically three state variables Ss,Xs and Xbh. The rate equations look like this:</p>
  
  <blockquote>
    <p>r1=µH(Ss/(Ks+Ss))(So/(Koh+So))Xbh+Kh(
    (Xs⁄Xbh)/(Xs⁄Xbh+Kx))(So/(Koh+So))Xbh</p>
    
    <p>r2=(1-fp)bH*Xbh-Kh( (Xs⁄Xbh)/(Xs⁄Xbh+Kx))(So/(Koh+So))Xbh </p>
    
    <p>r3=µH(Ss/(Ks+Ss))(So/(Koh+So))Xbh-bH*Xbh</p>
  </blockquote>
  
  <p>And the main differential equations derived from mole balance for CSTR are:</p>
  
  <blockquote>
    <p>dSs/dt = Q(Ss_in-Ss)+r1*V</p>
    
    <p>dXs/dt= Q(Xs_in-Xs)+r2*V</p>
    
    <p>dXbh/dt= Q(Xbh_in-Xbh)+r2*V</p>
  </blockquote>
  
  <p>Here is my code till now:</p>
  
  <pre><code>import numpy as np
  from scipy.optimize import fsolve
  parameter=dict()
  parameter['u_h']=6.0
  parameter['k_oh']=0.20
  parameter['k_s']=20.0
  parameter['k_h']=3.0
  parameter['k_x']=0.03
  parameter['Y_h']=0.67
  parameter['f_p']=0.08
  parameter['b_h']=0.62 
  
  Bulk_DO=2.0 #mg/L
  
  #influent components:
  infcomp=[56.53,182.9,16.625] #mgCOD/l
  
  Q=684000 #L/hr
  V=1040000 #l
  
  
  
  def steady(z,*args):
      Ss=z[0]
      Xs=z[1]
      Xbh=z[2]
      def monod(My_S,My_K):
          return My_S/(My_S+My_K)
  
      #Conversion rates
      #Conversion of Ss
      r1=((-1/parameter['Y_h'])*parameter['u_h']*monod(Ss,parameter['k_s'])\
          +parameter['k_h']*monod(Xs/Xbh,parameter['k_x'])*monod(Bulk_DO,parameter['k_oh']))\
          *Xbh*monod(Bulk_DO,parameter['k_oh'])
  
      #Conversion of Xs
      r2=((1-parameter['f_p'])*parameter['b_h']-parameter['k_h']*monod(Xs/Xbh,parameter['k_x']))*Xbh
  
      #Conversion of Xbh
      r3=(parameter['u_h']*monod(Ss,parameter['k_s'])*monod(Bulk_DO,parameter['k_oh'])-parameter['b_h'])*Xbh
  
      f=np.zeros(3)
      f[0]=Q*(infcomp[0]-Ss)+r1*V
      f[1]=Q*(infcomp[1]-Xs)+r2*V
      f[2]=Q*(infcomp[2]-Xbh)+r3*V
      return f
  initial_guess=(0.1,0.1,0.1)
  soln=fsolve(steady,initial_guess,args=parameter)
  print (soln)
  </code></pre>
  
  <p>How can I plot steady condition like this?
  <a href="https://i.stack.imgur.com/e5or9.jpg" rel="nofollow noreferrer">steady state plot</a>
  The solution is also not what I want since the equations implies reduction in Ss and Xs and increase of Xbh values with time. Also one solution has negative value which is practically impossible.
  Any suggestions would be highly appreciated. Thanks in advance !!</p>