code: def _z2deriv(self,R,z,phi=0.,t=0.):
          
          Raz2= (R+self.a)**2+z**2
          m= 4.*R*self.a/Raz2
          # Explicitly swapped in zforce here, so the z/z can be cancelled 
          # and z=0 is handled properly
          return -4.*(3.*z**2/Raz2-1.
                      +4.*((1.+m)/(1.-m)-special.ellipk(m)/special.ellipe(m))\
                          *self.a*R*z**2/Raz2**2/m)\
                      *self.a/(1.-m)*((R+self.a)**2+z**2)**-1.5*special.ellipe(m)
docs: NAME:
             _z2deriv
          PURPOSE:
             evaluate the second vertical derivative for this potential
          INPUT:
             R - Galactocentric cylindrical radius
             z - vertical height
             phi - azimuth
             t- time
          OUTPUT:
             the second vertical derivative
          HISTORY:
             2018-08-04 - Written - Bovy (UofT)