title: How to optimize a system of differential equations with curve_fit and solve_ivp?
body: <p>I want to optimize parameters in a system of differential equations. For a minimal example, I have the differential equations: <br /> <code>dx1/dt=-k*x1</code> and <br /> <code>dx2/dt=k*x1</code>.</p>
  <p>Therefore, I use <code>solve_ivp</code> to integrate the ODE's and <code>curve_fit</code> to fit the parameters.</p>
  <pre><code>t_end=5
  t_eval=np.linspace(0,t_end,t_end+1) #define timespan of 
  
  def fun(t,x,k): #System of differential equations with parameter k to optimize
         return [-k*x[0],k*x[0]]
  
  def fun2(t,k): #Solving differential equations and return x2(t)
         res=solve_ivp(fun, [0,t_end], [1,0],t_eval=t_eval,args=(k,))
         return res.y[1]
  
  z=(1-np.exp(-0.2*t_eval))*np.ones(len(t_eval)) #creating an array with the analytic solution for the above statet differential equations with k=0.2
  popt, pcov = curve_fit(fun2, t_eval, z)
  </code></pre>
  <p>This works, but if I increase <code>t_end</code> or make the differential equation more complicated, it always gets the same errors:</p>
  <blockquote>
  <p>ValueError: operands could not be broadcasted together with shapes (xxx,) (xxx,)</p>
  </blockquote>
  <p>e.g. with t_end=1200 i get the error</p>
  <blockquote>
  <p>ValueError: operands could not be broadcasted together with shapes (857,) (1201,)</p>
  </blockquote>
  <p>or, with the slightly more complex system of ODE's and <code>t_end = 10</code>:</p>
  <pre><code>def fun(t,x,k1,k2): #System of differential equations
     return [-k1*x[0],k1*x[0]-k2*x[1],k2*x[1]]
  
  def fun2(t,k1,k2): #Solving differential equations and return x3(t)
     res=solve_ivp(fun, [0,t_end], [1,0,0],t_eval=t_eval,args=(k1,k2))
     return res.y[2]
  </code></pre>
  <p>I get the below error:</p>
  <blockquote>
  <p>ValueError: operands could not be broadcasted together with shapes (8,) (11,)</p>
  </blockquote>
  <p>I already figured out two options in <code>solve_ivp</code> which let me use higher <code>t_end</code>'s and more complex systems of ODE's: Choosing <code>method='DOP853'</code> or <code>method='LSODA'</code> with <code>min_step=1e-3</code> or <code>1e-2</code>.</p>
  <p>But in the end, it's not enough for the System I want to use (4-5 differential equations with 7-8 parameters to fit and t_end of 6000+).</p>
  <p>So, Is there any general problem with my approach or any trick I missed?</p>