'This function does compute series for multivariate functions,
but the expansion is always in terms of *one* variable.
Examples
>>> from sympy import atan2
>>> from sympy.abc import x, y
>>> atan2(x, y).series(x, n=2)
atan2(0, y) + x/y + O(x**2)
>>> atan2(x, y).series(y, n=2)
-y/x + atan2(x, 0) + O(y**2)
This function also computes asymptotic expansions, if necessary
and possible:
>>> from sympy import loggamma
>>> loggamma(1/x)._eval_nseries(x,0,None)
-1/x - log(x)/x + log(x)/2 + O(1)'

def _eval_nseries(self, x, n, logx):

    from sympy import Order
    from sympy.sets.sets import FiniteSet
    args = self.args
    args0 = [t.limit(x, 0) for t in args]
    if any(((t.is_finite is False) for t in args0)):
        from sympy import oo, zoo, nan
        a = [t.compute_leading_term(x, logx=logx) for t in args]
        a0 = [t.limit(x, 0) for t in a]
        if any([t.has(oo, (- oo), zoo, nan) for t in a0]):
            return self._eval_aseries(n, args0, x, logx)
        a = [t._eval_nseries(x, n, logx) for t in args]
        z = [(r - r0) for (r, r0) in zip(a, a0)]
        p = [Dummy() for t in z]
        q = []
        v = None
        for (ai, zi, pi) in zip(a0, z, p):
            if zi.has(x):
                if (v is not None):
                    raise NotImplementedError
                q.append((ai + pi))
                v = pi
            else:
                q.append(ai)
        e1 = self.func(*q)
        if (v is None):
            return e1
        s = e1._eval_nseries(v, n, logx)
        o = s.getO()
        s = s.removeO()
        s = (s.subs(v, zi).expand() + Order(o.expr.subs(v, zi), x))
        return s
    if ((self.func.nargs is S.Naturals0) or ((self.func.nargs == FiniteSet(1)) and args0[0]) or any(((c > 1) for c in self.func.nargs))):
        e = self
        e1 = e.expand()
        if (e == e1):
            term = e.subs(x, S.Zero)
            if ((term.is_finite is False) or (term is S.NaN)):
                raise PoleError(('Cannot    expand    %s    around    0' % self))
            series = term
            fact = S.One
            _x = Dummy('x')
            e = e.subs(x, _x)
            for i in range((n - 1)):
                i += 1
                fact *= Rational(i)
                e = e.diff(_x)
                subs = e.subs(_x, S.Zero)
                if (subs is S.NaN):
                    subs = e.limit(_x, S.Zero)
                if (subs.is_finite is False):
                    raise PoleError(('Cannot    expand    %s    around    0' % self))
                term = ((subs * (x ** i)) / fact)
                term = term.expand()
                series += term
            return (series + Order((x ** n), x))
        return e1.nseries(x, n=n, logx=logx)
    arg = self.args[0]
    l = []
    g = None
    nterms = (n + 2)
    cf = Order(arg.as_leading_term(x), x).getn()
    if (cf != 0):
        nterms = int((nterms / cf))
    for i in range(nterms):
        g = self.taylor_term(i, arg, g)
        g = g.nseries(x, n=n, logx=logx)
        l.append(g)
    return (Add(*l) + Order((x ** n), x))
