title: Python: How to solve an ordinary differential equation with integral term in it
body: <p>I dont know if this question has been asked before in SO, I will go ahead and post it here, I am attempting to solve a simple system with a PID controller, my system of differential equations are given below. I am basically attempting to code very basic PID algorithm. The structure of my control u depends on both derivative and integral of error term. I dont have any problem with the derivative term, it is the integral term that is creating problem in my code. The problem crops up when I assign s=0 in the beginning 
  and use it in my function as described in my code below. Is there a way to bypass it? I tried assigning s and told as global variables, but it didnt solve my problem. In a nutshell what I am doing is- I am adding state x1 every time and multiplying by dt(which is denoted by t-told).</p>
  
  <p><a href="https://i.stack.imgur.com/eHDZf.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/eHDZf.png" alt="enter image description here"></a></p>
  
  <p>Kindly help me iron out this issue, PFA my code attached below.</p>
  
  <pre><code>import numpy as np
  from scipy.integrate import ode
  import matplotlib.pyplot as plt
  plt.style.use('bmh')
  
  t0=0
  y0=[0.1,0.2]
  kp,kd,ki=2,0.5,0.8
  s,told=0,0
  def pid(t,Y):
      x1,x2=Y[0],Y[1]
      e=x1-1
      de=x2
      s=(x1+s)
      integral=s*(t-told)
      told=t
      #ie=
      u=kp*e+kd*de+ki*integral
      x1dot=x2
      x2dot=u-5*x1-2*x2
      return[x1dot,x2dot]
  
  solver=ode(pid).set_integrator('dopri5',rtol=1e-6,method='bdf',nsteps=1e5,max_step=1e-3)
  solver.set_initial_value(y0,t0)
  t1=10
  dt=5e-3
  sol = [ [yy] for yy in y0 ]
  t=[t0]
  while solver.successful() and solver.t&lt;t1:
      solver.integrate(solver.t+dt)
      for k in range(2): sol[k].append(solver.y[k]);
      t.append(solver.t)
      print(len(sol[0]))
      print(len(t))
  x1=np.array(sol[0])
  x2=np.array(sol[1])
  e=x1-1
  de=x2
  u=kp*e+kd*de
  for k in range(2):
      if k==0:
          plt.subplot(2,1,k+1)
          plt.plot(t,sol[k],label='x1')
          plt.plot(t,sol[k+1],label='x2')
          plt.legend(loc='lower right')
      else:
          plt.subplot(2,1,k+1)
          plt.plot(t,u)
  plt.show()
  </code></pre>