title: Does the order of the equations in a coupled OIDENT solver matter?
body: <p>My code runs, however, if i change the order of the equations defined in my formula my graphs changes as well. Can somebody tell me why this is? Now i do not know what the right graph is from this system. </p>
  
  <pre><code>def myFunction(r,t):
  g = 9.81
  L_L = 20 #draught
  L_r = 20 #draught
  L_d = 4 #ukc
  u_s = 0.08
  w_d = 60 #width vessel
  rho = 1030
  b_R = 1.0
  b_L = 1.0
  b_D = 1.0
  A_s = 340*L_L
  M_s = 40000*10^3 
  
  w_L = r[0]
  w_r = r[1]
  u_1 = r[2]
  u_2 = r[3]
  u_3 = r[4]
  u_4 = r[5]
  u_d = r[6]
  p_1 = r[7]
  p_2 = r[8]
  p_3 = r[9]
  p_4 = r[10]
  deltap_L = r[11]
  deltap_R = r[12]
  
  u_1 = (L_L*u_s + w_L*u_2)/w_L
  u_2 = (w_d*u_d)/w_L
  u_d = (w_r*u_3)/w_d
  u_3 = (w_L*u_2)/w_r
  u_4 = (L_r*u_s + w_r*u_3)/w_r
  dwLdt = u_s
  dwrdt = - u_s  
  du1dt =  - g*((p_2-p_1)/L_L + deltap_L/L_L) - b_L*u_1 
  du2dt =  - g*((p_2-p_1)/L_L - deltap_L/L_L) - b_L*u_2
  du3dt =  - g*((p_4-p_3)/L_r + deltap_R/L_r) - b_R*u_3
  du4dt =  - g*((p_4-p_3)/L_r - deltap_R/L_r) - b_R*u_4
  duddt =  - g*((p_3-p_2)/L_d) - b_D *u_d
  dp1dt = - u_1
  dp2dt = - u_1
  dp3dt = + u_4
  dp4dt = + u_4
  ddeltap_Ldt = -u_1
  ddeltap_Rdt = u_4 
  #M_s = (g*rho*A_s*((p_1+p_2)/2 + deltap_L/6 - (p_3+p_4)/2 - deltap_R/6))/u_s
  
  
  return (dwLdt, dwrdt, du1dt, du2dt, du3dt, du4dt, duddt, dp1dt, dp2dt, dp3dt, dp4dt, ddeltap_Ldt, ddeltap_Rdt)
  r0 = [10,10,0,0,0,0,0,0,0,0,0,0,0]
  t = np.linspace(0,20,10000)
  r = odeint(myFunction, r0, t)
  
  w_L = r[:,0]
  w_r = r[:,1]
  u1 = r[:,2]
  u2 = r[:,3]
  u3 = r[:,4]
  u4 = r[:,5]
  ud = r[:,6]
  p_1 = r[:,7]
  p_2 = r[:,8]
  p_3 = r[:,9]
  p_4 = r[:,10]
  deltap_L = r[:,11]
  deltap_R = r[:,12]
  
  plt.figure()
  plt.plot(t, u1)
  plt.plot(t, u2)
  plt.plot(t, u3)
  plt.plot(t, u4)
  plt.plot(t, ud)
  plt.legend(('$u_1$', '$u_2$', '$u_3$', '$u_4$', '$u_d$'))
  plt.xlabel('Time ($s$)')
  plt.ylabel('Fluid velocity ($m/s$)')
  </code></pre>
  
  <p><a href="https://i.stack.imgur.com/OtC5e.png" rel="nofollow noreferrer">graph first code</a></p>
  
  <p>if the equations in the code are in this order:</p>
  
  <pre><code>    dwLdt = u_s
      dwrdt = - u_s  
      du1dt =  - g*((p_2-p_1)/L_L + deltap_L/L_L) - b_L*u_1 
      du2dt =  - g*((p_2-p_1)/L_L - deltap_L/L_L) - b_L*u_2
      du3dt =  - g*((p_4-p_3)/L_r + deltap_R/L_r) - b_R*u_3
      du4dt =  - g*((p_4-p_3)/L_r - deltap_R/L_r) - b_R*u_4
      duddt =  - g*((p_3-p_2)/L_d) - b_D *u_d
      u_1 = (L_L*u_s + w_L*u_2)/w_L
      u_2 = (w_d*u_d)/w_L
      u_d = (w_r*u_3)/w_d
      u_3 = (w_L*u_2)/w_r
      u_4 = (L_r*u_s + w_r*u_3)/w_r
      dp1dt = - u_1
      dp2dt = - u_1
      dp3dt = + u_4
      dp4dt = + u_4
      ddeltap_Ldt = -u_1
      ddeltap_Rdt = u_4 
      #M_s = (g*rho*A_s*((p_1+p_2)/2 + deltap_L/6 - (p_3+p_4)/2 - deltap_R/6))/u_s
  </code></pre>
  
  <p>the graph is like this:</p>
  
  <p><a href="https://i.stack.imgur.com/U4eAI.png" rel="nofollow noreferrer">graph second mode</a></p>
  
  <p>Why is this the case?</p>
  
  <p>The code is now:</p>
  
  <pre><code>  def myFunction(r,t):
       g = 9.81
       L_L = 20 #draught
       L_r = 20 #draught
       L_d = 4 #ukc
       u_s = 0.08
       w_d = 60 #width vessel
       rho = 1025
       b_R = 1.0
       b_L = 1.0
       b_D = 1.0
       A_s = 340*L_L
       M_s = 40000*10^6
  
       w_L = r[0]
       w_r = r[1]
       u_1 = r[2]
       u_2 = r[3]
       u_3 = r[4]
       u_4 = r[5]
       u_d = r[6]
       p_1 = r[7]
       p_2 = r[8]
       p_3 = r[9]
       p_4 = r[10]
       deltap_L = r[11]
       deltap_R = r[12]
  
  
       u_1 = (L_L*u_s + w_L*u_2)/w_L
       u_2 = (w_d*u_d)/w_L
       u_3 = (w_L*u_2)/w_r
       u_4 = (L_r*u_s + w_r*u_3)/w_r
       u_d = (w_r*u_3)/w_d
  
       du1dt =  - g*((p_2-p_1)/L_L + deltap_L/L_L) - b_L*u_1 
       du2dt =  - g*((p_2-p_1)/L_L - deltap_L/L_L) - b_L*u_2
       du3dt =  - g*((p_4-p_3)/L_r - deltap_R/L_r) - b_R*u_3
       du4dt =  - g*((p_4-p_3)/L_r + deltap_R/L_r) - b_R*u_4
       duddt =  - g*((p_3-p_2)/L_d) - b_D *u_d
  
  
  
  
       dp1dt =  - u_1
       dp2dt =  - u_1
       dp3dt = + u_4
       dp4dt = + u_4
       ddeltap_Ldt =  - u_1
       ddeltap_Rdt =  u_4 
  
       dwLdt = u_s
       dwrdt = - u_s 
  
       deltap_R = 6*(((p_1+p_2)/2) + (deltap_L/6) - ((p_3+p_4)/2)) 
  
       return (dwLdt, dwrdt, du1dt, du2dt, du3dt, du4dt, duddt, dp1dt,     dp2dt, dp3dt, dp4dt, ddeltap_Ldt, ddeltap_Rdt)
  
  r0 = [10,10,0,0,0,0,0,0,0,0,0,0,0]
  t = np.linspace(0,125,100000)
  r = odeint(myFunction, r0, t)
  
  wL = r[:,0]
  wr = r[:,1]
  u1 = r[:,2]
  u2 = r[:,3]
  u3 = r[:,4]
  u4 = r[:,5]
  ud = r[:,6]
  p1 = r[:,7]
  p2 = r[:,8]
  p3 = r[:,9]
  p4 = r[:,10]
  deltapL = r[:,11]
  deltapR = r[:,12]
  </code></pre>