title: "RuntimeError: The array return by func must be one-dimensional, but got ndim=2."
body: <p>I want to solve this 3-3 differential equation system using odeint, and I get this error: </p>
  
  <pre><code>Traceback (most recent call last):
    File "/home/spyros/Documents/Spyros/Diplwmatiki/ROS/Serial connection/Dynamic_problem.py", line 81, in &lt;module&gt;
      u = odeint(subsystem2,u0,t)
    File "/usr/local/lib/python2.7/dist-packages/scipy/integrate/odepack.py", line 233, in odeint
      int(bool(tfirst)))
  RuntimeError: The array return by func must be one-dimensional, but got ndim=2.
  </code></pre>
  
  <p>I read these similar topics: 
  <a href="https://stackoverflow.com/questions/51808922/how-to-solve-a-system-of-differential-equations-using-scipy-odeint">How to solve a system of differential equations using scipy.odeint</a>
  <a href="https://stackoverflow.com/questions/51842875/error-with-odeint-runtimeerror-the-array-return-by-func-must-be-one-dimensiona">error with odeint (RuntimeError: The array return by func must be one-dimensional, but got ndim=2.)</a>
  <a href="https://stackoverflow.com/questions/51580980/scipy-optimize-error">Scipy optimize error</a>
  but I didn't find how the solutions can be implement in my problem.
  How is it possible that the returning array has dimension = 2. I have 3 variables in the function's input, 3 differential equations and the same equations in the returning array. How can I solve this error? Thank you in advance.</p>
  
  <pre><code>import numpy as np 
  from scipy.integrate import odeint 
  import matplotlib.pyplot as plt
  
  added_mass_x = 0.03 # kg 
  added_mass_y = 0.04 
  mb = 0.3 # kg 
  m1 = mb-added_mass_x 
  m2 = mb-added_mass_y 
  l1 = 0.07 # m 
  l2 = 0.05 # m 
  J = 0.00050797 # kgm^2 
  Sa = 0.0110 # m^2 
  Cd = 2.44 
  Cl = 3.41 
  Kd = 0.000655 # kgm^2 
  r = 1000 # kg/m^3
  
  c1 = 0.5*r*Sa*Cd 
  c2 = 0.5*r*Sa*Cl 
  c3 = 0.5*mb*(l1**2) 
  c4 = Kd/J 
  c5 = (1/(2*J))*(l1**2)*mb*l2 
  c6 = (1/(3*J))*(l1**3)*mb
  
  theta_0 = 10*(np.pi/180) # rad 
  theta_A = 20*(np.pi/180) # rad 
  f = 2 # Hz
  
  t = np.linspace(0,100,8000) # s 
  omega = 2*np.pi*f # rad/s 
  theta = theta_0 + theta_A*np.sin(omega*t) # rad 
  theta_deg = (theta*180)/np.pi # deg 
  thetadotdot = -(omega**2)*theta_A*np.sin(omega*t) # rad/s^2
  
  def subsystem2(u,t):
      vcx = u[0]
      vcy = u[1]
      wz = u[2]
  
      vcxdot = (m2/m1)*vcy*wz-(c1/m1)*vcx*np.sqrt((vcx**2)+(vcy**2))+(c2/m1)*vcy*np.sqrt((vcx**2)+(vcy**2))*np.arctan2(vcy,vcx)-(c3/m1)*thetadotdot*np.sin(theta)
      vcydot = -(m1/m2)*vcx*wz-(c1/m2)*vcy*np.sqrt((vcx**2)+(vcy**2))-(c2/m2)*vcx*np.sqrt((vcx**2)+(vcy**2))*np.arctan2(vcy,vcx)+(c3/m2)*thetadotdot*np.cos(theta)
      wzdot = ((m1-m2)/J)*vcx*vcy-c4*wz*wz*np.sign(wz)-c5*thetadotdot*np.cos(theta)-c6*thetadotdot
  
      return [vcxdot,vcydot,wzdot]
  
  # arxikes synthikes 
  vcx_0 = 0 
  vcy_0 = 0 
  wz_0 = 0 
  u0 = [vcx_0,vcy_0,wz_0]
  
  # epilysi systimatos diaforikwn 
  u = odeint(subsystem2,u0,t,tfirst=False)
  
  vcx = u[:,0] 
  vcy = u[:,1] 
  wz = u[:,2]
  </code></pre>