code: def _R2deriv(self,R,z,phi=0.,t=0.):
          
          Raz2= (R+self.a)**2+z**2
          Raz= nu.sqrt(Raz2)
          m= 4.*R*self.a/Raz2
          R2ma2mz2o4aR1m= (R**2-self.a2-z**2)/4./self.a/R/(1.-m)
          return (2*R**2+self.a2+3*R*self.a+z**2)/R/Raz2*self._Rforce(R,z)\
              +2.*self.a/R/Raz*(m*(R**2+self.a2+z**2)/4./(1.-m)/self.a/R**2\
                                    *special.ellipe(m)\
                +(R2ma2mz2o4aR1m/(1.-m)*special.ellipe(m)
                  +0.5*R2ma2mz2o4aR1m*(special.ellipe(m)-special.ellipk(m))
                  +0.5*(special.ellipe(m)/(1.-m)-special.ellipk(m))/m)\
                                    *4*self.a*(self.a2+z**2-R**2)/Raz2**2)
docs: NAME:
             _Rderiv
          PURPOSE:
             evaluate the second radial derivative for this potential
          INPUT:
             R - Galactocentric cylindrical radius
             z - vertical height
             phi - azimuth
             t - time
          OUTPUT:
             the second radial derivative
          HISTORY:
             2018-08-04 - Written - Bovy (UofT)