code: def _Rzderiv(self,R,z,phi=0.,t=0.):
          
          l,n    = bovy_coords.Rz_to_lambdanu    (R,z,ac=self._ac,Delta=self._Delta)
          jac    = bovy_coords.Rz_to_lambdanu_jac(R,z,            Delta=self._Delta)
          hess   = bovy_coords.Rz_to_lambdanu_hess(R,z,            Delta=self._Delta)
          dldR = jac[0,0]
          dndR = jac[1,0]
          dldz = jac[0,1]
          dndz = jac[1,1]
          d2ldRdz = hess[0,0,1]
          d2ndRdz = hess[1,0,1]
          return d2ldRdz              * self._lderiv(l,n)  + \
                 d2ndRdz              * self._nderiv(l,n)  + \
                 dldR*dldz            * self._l2deriv(l,n) + \
                 dndR*dndz            * self._n2deriv(l,n) + \
                 (dldR*dndz+dldz*dndR)* self._lnderiv(l,n)
docs: NAME:
              _Rzderiv
          PURPOSE:
              evaluate the mixed R,z derivative for this potential
          INPUT:
              R - Galactocentric cylindrical radius
              z - vertical height
              phi - azimuth
              t- time
          OUTPUT:
              d2phi/dR/dz
          HISTORY:
              2015-02-13 - Written - Trick (MPIA)