title: scipy optimization solve non linear equations with array?
body: <p>I have to solve the follow sistem of equation :</p>
  
  <pre><code>f1 = -k1 + self.f(
      self.time[i] + (0.5 + np.sqrt(3)/6) * self.dt,
      self.u[i] + 0.25 * self.dt * k1 + (0.25 + np.sqrt(3)/6) * self.dt * k2) 
  f2 = -k2 + self.f(
      self.time[i] + (0.5 - np.sqrt(3)/6) * self.dt,
      self.u[i] + (0.25 - np.sqrt(3)/6) * self.dt * k1 + 0.25 * self.dt * k2)
  </code></pre>
  
  <p>Searching in the previous post I found this solutions :</p>
  
  <pre><code>def equations(variable):
      k1,k2 = variable
      f1 = -k1 + self.f(
          self.time[i] + (0.5 + np.sqrt(3)/6) * self.dt,
          self.u[i] + 0.25 * self.dt * k1 + (0.25 + np.sqrt(3)/6) * self.dt * k2)
      f2 = -k2 + self.f(
          self.time[i] + (0.5 - np.sqrt(3)/6) * self.dt,
          self.u[i] + (0.25 - np.sqrt(3)/6) * self.dt * k1 + 0.25 * self.dt * k2)
      return np.array([f1,f2]).ravel()     #.reshape(2,)  
  
  k1, k2 = fsolve(equations, (self.u[i],self.u[i]))
  </code></pre>
  
  <p>But does't works for me, because my <code>self.u[i]</code> in on <code>len(self.u[i]) == 2</code> 
  how can I deal in this situations??</p>
  
  <p>I got this error:</p>
  
  <pre><code>    k1,k2 = variable
  ValueError: too many values to unpack (expected 2)
  </code></pre>
  
  <p>but I don't know how to fix it , if i simply pass 2 float <code>(2.,2.)</code> I got  this :</p>
  
  <pre><code>TypeError: fsolve: there is a mismatch between the input and output shape of the 'func' argument 'equations'.
  Shape should be (2,) but it is (4). 
  </code></pre>
  
  <p>and if i left the <code>.ravel()</code> on the return array i got the same error but :</p>
  
  <pre><code>TypeError: fsolve: there is a mismatch between the input and output shape of the 'func' argument 'equations'.
  Shape should be (2,) but it is (2,2). 
  </code></pre>
  
  <p>please may somebody help me to come up with this ?</p>
  
  <pre><code>class ImpRK4 :
      def __init__(self, fun , t0, tf, dt , y0):
          self.func = fun
          self.t0=t0
          self.tf=tf
          self.dt=dt
          self.u0=y0
          self.n = round((tf-t0)/dt)
          self.time  = np.linspace(self.t0, self.tf, self.n+1 )
          self.u     = np.array([self.u0  for i in range(self.n+1) ])
  
      def f(self,ti,ui):
           return  np.array([functions(ti,ui) for functions in self.func])     
  
      def solve(self): 
         for i in range(len(self.time)-1):
              def equations(variable):
                  k1,k2 = variable
                  f1 = -k1 + self.f(self.time[i]+ (0.5+np.sqrt(3)/6)* self.dt , self.u[i]+0.25*self.dt* k1+ (0.25+ np.sqrt(3)/6)*self.dt*k2) 
                  f2 = -k2 + self.f(self.time[i]+ (0.5-np.sqrt(3)/6)* self.dt , self.u[i]+(0.25-np.sqrt(3)/6)*self.dt *k1 + 0.25*self.dt* k2)
                  return np.array([f1,f2]).ravel() #.reshape(2,)  
  
              k1 , k2 = fsolve(equations,(2,2)) #(self.u[i],self.u[i]))
              self.u[i+1] = self.u[i] + self.dt/2* (k1 + k2)
  
         plt.plot(self.time,self.u)
         plt.show()    
  
  def main():
      func00 = lambda t,u : -10*(t-1)*u[0]
  
      func01 = lambda t,u : u[1] 
      func02 = lambda t,u : (1-u[0]**2)*u[1] - u[0]
  
      func0x = np.array([func00])
      func0 = np.array([func01,func02])
  
      y01 = np.array([1.,1.])
      t0 = 0. 
      tf = 2.      
      u0 = y01   
      dt = 0.008 
  
      y0  = np.array([np.exp(-5)])
      diffeq = ImpRK4(func0,t0,tf,dt,y01)    
  
      #diffeq.solve()
      #diffeq = ImpRK4(func0x,t0,tf,dt,y0) ## with single equations works
      diffeq.solve()
  
  if __name__ == '__main__': 
      main()
  </code></pre>
  
  <p>Ok here the Minimal, Complete, and Verifiable example</p>