Compute a generalized series expansion

MuPAD® notebooks will be removed in a future release. Use MATLAB® live scripts instead.

MATLAB live scripts support most MuPAD functionality, though there are some differences. For more information, see Convert MuPAD Notebooks to MATLAB Live Scripts.


series(f, x, <order>, <Left | Right | Real | Undirected>, <NoWarning>, <UseGseries>)
series(f, x = x0, <order>, <Left | Right | Real | Undirected>, <options>)


series(f, x = x0) computes the first terms of a series expansion of f with respect to the variable x around the point x0.

series tries to compute either the Taylor series, the Laurent series, the Puiseux series, or a generalized series expansion of f around x = x0. See Series::gseries for details on generalized series expansions.

The mathematical type of the series returned by series can be queried using the type expression Type::Series.

If series cannot compute a series expansion of f, a symbolic function call is returned. This is an expression of type "series". Cf. Example 11.

Mathematically, the expansion computed by series is valid in some neighborhood of the expansion point in the complex plane. Usually, this is an open disc centered at x0. However, if the expansion point is a branch point, then the returned expansion may not approximate the function f for values of x close to the branch cut. Cf. Example 12.

Using the options Left or Right, one can compute directed expansions that are valid along the real axis. With the option Real, a two-sided expansion along the real axis is computed. See Example 5 and Example 6.

If x0 is infinity or -infinity, then a directed series expansion along the real axis from the left to the positive real infinity or from the right to the negative real infinity, respectively, is computed. If x0 is complexInfinity and dir is not specified or Undirected, then an undirected series expansion around the complex infinity, i.e., the north pole of the Riemann sphere, is computed. Specifying x0= infinity is equivalent to x0= complexInfinity and dir = Left. Similarly, x0= -infinity is equivalent to x0= complexInfinity and dir = Right. Cf. Example 7.

Such a series expansion is computed as follows: The series variable x in f is replaced by (or for x0= -infinity). Then, a series expansion of f around u = 0 is computed. Finally, (or , respectively) is substituted in the result.

Mathematically, the result of such a series expansion is a series in . However, it may happen that the coefficients of the returned series depend on the series variable. See the corresponding paragraph below.

The number of requested terms for the expansion is the argument order if specified. Otherwise, the value of the environment variable ORDER is used. One can change the default value 6 by assigning a new value to ORDER.

The number of terms is counted from the lowest degree term on for finite expansion points, and from the highest degree term on for expansions around infinity, i.e., “order” has to be regarded as a “relative truncation order”.

series implements a limited amount of precision management to circumvent cancellation. If the number of terms of the computed expansion is less than order, a second series computation with a higher value of order is tried automatically, and the result of the latter is returned.


Nevertheless, the actual number of terms in the resulting series expansion may differ from the requested number of terms. See Example 13 and Example 15.

Taylor/Laurent/Puiseux expansions (all of domain type Series::Puiseux) can be restricted easily to an absolute order term by adding an appropriate O term. Cf. Example 14.

Expansions of symbolic integrals can be computed. Cf. Example 16.

If f is an expression of type RootOf, then series returns the set of all non-zero series solutions of the corresponding algebraic equation. Cf. Example 9.

If order has the value infinity, then the system tries to convert the first argument into a formal infinite series, i.e., it computes a general formula for the n-th coefficient in the Taylor expansion of f. The result is an inactive symbolic sum or a polynomial expression. Cf. Example 10.

If series returns a series expansion of domain type Series::Puiseux, it may happen that the “coefficients” of the returned series depend on the series variable. In this case, the expansion is not a proper Puiseux series in the mathematical sense. See Example 7 and Example 8. However, if the series variable is x and the expansion point is x0, then the following is valid for each coefficient function c(x) and every positive ε: c(x) (x - x0)ε converges to zero and is unbounded when x approaches x0. Similarly, if the expansion point is infinity, then, for every positive ε, converges to zero and c(x) xε is unbounded when x approaches infinity.

The function returns a domain object that can be manipulated by the standard arithmetical operations. Moreover, the following methods are available: ldegree returns the exponent of the leading term; Series::Puiseux::order returns the exponent of the error term; expr converts to an arithmetical expression, removing the error term; coeff(s, n) returns the coefficient of the term of s with exponent n; lcoeff returns the leading coefficient; revert computes the inverse with respect to composition; diff and int differentiate and integrate a series expansion, respectively; map applies a function to all coefficients. See the help pages for Series::Puiseux and Series::gseries for further details.


series works on a symbolic level and should not be called with arguments containing floating point arguments.

Environment Interactions

The function is sensitive to the environment variable ORDER, which determines the default number of terms in series computations.


Example 1

We compute a series expansion of sin(x) around x = 0. The result is a Taylor series:

s := series(sin(x), x)

Syntactically, the result is an object of domain type Series::Puiseux:


The mathematical type of the series expansion can be queried using the type expression Type::Series:

testtype(s, Type::Series(Taylor))

Various system functions are overloaded to operate on series objects. E.g., the function coeff can be used to extract the coefficients of a series expansion:

coeff(s, 5)

The standard arithmetical operators can be used to add or multiply series expansions:

s + 2*s, s*s

delete s:

Example 2

This example computes the composition of s by itself, i.e. the series expansion of sin(sin(x)).

s := series(sin(x), x): s @ s = series(sin(sin(x)), x)

delete s:

Example 3

We compute the series expansion of the tangent function around the origin in two ways:

series(sin(x), x) / series(cos(x), x) = series(tan(x), x)


Example 4

We compute a Laurent expansion around the point 1:

s := series(1/(x^2 - 1), x = 1)

testtype(s, Type::Series(Taylor)),
testtype(s, Type::Series(Laurent))

Example 5

Without an optional argument or with the option Undirected, the sign function is not expanded:

series(x*sign(x^2 + x), x) =
series(x*sign(x^2 + x), x, Undirected)

Some simplification occurs if one requests an expansion that is valid along the real axis only:

series(x*sign(x^2 + x), x, Real)

The sign vanishes from the result if one requests a one-sided expansion along the real axis:

series(x*sign(x^2 + x), x, Right),
series(x*sign(x^2 + x), x, Left)

Example 6

In MuPAD®, the heaviside function is defined only on the real axis. Thus an undirected expansion in the complex plane does not make sense:

series(x*heaviside(x + 1), x)
Warning: Unable to find an undirected series expansion. Try the 'Left', 'Right', or 'Real' option. [Series::main]

After specifying corresponding options, the system computes an expansion along the real axis:

series(x*heaviside(x + 1), x, Real),
series(x*heaviside(x + 1), x, Right)

At the point I in the complex plane, the function heaviside is not defined, and neither is a series expansion:

series(heaviside(x), x = I, Real)

Example 7

We compute series expansions around infinity:

s1 := series((x + 1)/(x - 1), x = complexInfinity)

s2 := series(psi(x), x = infinity)

domtype(s1), domtype(s2)

Although both expansions are of domain type Series::Puiseux, s2 is not a Puiseux series in the mathematical sense, since the first term contains a logarithm, which has an essential singularity at infinity:

testtype(s1, Type::Series(Puiseux)),
testtype(s2, Type::Series(Puiseux))


The following expansion is of domain type Series::gseries:

s3 := series(exp(x)/(1 - x), x = infinity, 4)


delete s1, s2, s3:

Example 8

Oscillating but bounded functions may appear in the “coefficients” of a series expansion as well:

s := series(sin(x + 1/x), x = infinity)

domtype(s), testtype(s, Type::Series(Puiseux))

coeff(s, -1)

Example 9

The algebraic equation y5 - y - x = 0 cannot be resolved in terms of radicals:

solve(y^5 - y - x, y)

However, series can compute all series solutions of this equation around x = 0:

series(%, x = 0)

It may happen that the series solutions themselves are expressed in terms of RootOfs:

series(RootOf(y^5 -(x + 2*x^2)*y^3 - x^3*y^2
              + (x^3 + x^4)*y + x^4 + x^5, y), x)

The coefficients of the algebraic equation are allowed to be transcendental. They are internally converted into Puiseux series by series:

series(RootOf(y^3 - y - exp(x - 1) + 1, y), x = 1, 4)

An error occurs if some coefficient cannot be expanded into a Puiseux series:

series(RootOf(y^3 - y - exp(x), y), x = infinity)
Error: Unable to expand coefficients of 'RootOf(y^3 - y - exp(1/x), y)' into a series. [Series::algebraic]

Example 10

In this example, we compute a formula for the n-th coefficient an in the Taylor expansion of the function around zero, by specifying infinity as order. The result is a symbolic sum:

series(exp(-x), x, infinity)

If the input is a polynomial expression, then so is the output:

series(x^5 - 1, x = 1, infinity)

Example 11

No asymptotic expansion exists for the Bessel J function of unspecified index, and series returns a symbolic function call:

series(besselJ(k, x), x=infinity)

domtype(%), type(%)

Example 12

The branch cut of the logarithm and the square root is the negative real axis. For a series expansion on the branch cut, series uses the function signIm to return an expansion that is valid in an open disc around the expansion point:

series(ln(x), x = -1, 3)

series(sqrt(x), x = -1, 3)

The situation is more intricate when the expansion point is a branch point. The following expansion of the function arcsin(x + 1) is valid in an open disc around the branch point 0:

series(arcsin(x + 1), x, 4)

However, the expansion of f = ln(x + I*x^3) around the branch point 0 that is returned by series does not approximate f for values of x that are close to the negative real axis:

f := ln(x + I*x^3);
g := series(f, x, 4);

DIGITS := 20:
float(subs([f, expr(g)], x = -0.01 + 0.0000001*I));
delete DIGITS:

The situation is similar for algebraic branch points:

f := sqrt(x + I*x^3);
g := series(f, x, 4);

DIGITS := 20:
float(subs([f, expr(g)], x = -0.01 + 0.0000001*I));
delete DIGITS:

delete f, g:

Example 13

The first six terms, including zeroes, of the following two series expansions agree:

series(sin(tan(x)), x, 12);
series(tan(sin(x)), x, 12);

If we want to compute the series expansion of the difference sin(tan(x)) - tan(sin(x)), cancellation happens and produces too few terms in the result. series detects this automatically and performs a second series computation with increased precision:

series(sin(tan(x)) - tan(sin(x)), x, 6)

It may nevertheless happen that the result has too few terms; cf. Example 15.

If rational exponents occur in the series expansion, then it may even happen that the result has more than the number of terms requested by the third argument:

series(x^2*exp(x) + x*sqrt(sin(x)), x, 3)

Example 14

series's control of the order term is based on the concept of `relative order', counting the number of terms beginning with the lowest order that is present in the expansion. An `absolute order' control can be achieved by simply adding an appropriate order term to restrict a result returned by series:

series(exp(x) + x*sqrt(sin(x)), x, 7)

series(exp(x) + x*sqrt(sin(x)), x, 7) + O(x^4)

Note, however, that the series must have enough terms for the added order term to have any effect:

series(exp(x^2), x, 4)

series(exp(x^2), x, 4) + O(x^8)

Example 15

If the specified order for the expansion is too small to compute the reciprocal (due to cancellation), series returns a symbolic call:

series(exp(x), x, 4)

series(1/(exp(x) - 1 - x - x^2/2 - x^3/6), x, 2)

After increasing the order, an expansion is computed, but possibly with fewer terms:

series(1/(exp(x) - 1 - x - x^2/2 - x^3/6), x, 3);
series(1/(exp(x) - 1 - x - x^2/2 - x^3/6), x, 4)

Example 16

series and int support each other. On the one hand, series expansions can be integrated:

int(series(1/(2 - x), x), x = 0..1)

On the other hand, series knows how to handle symbolic integrals:

int(x^x, x)

series(%, x = 0, 3)

int(exp(-x*sin(t)), t = 0.. x)

series(%, x = 0)

int(cos((x*t^2 + x^2*t))^(1/3), t = 0..2)

series(%, x)

Example 17

Users can extend the power of series by implementing series attributes (slots) for their own special mathematical functions.

We illustrate how to write such a series attribute, using the case of the exponential function. (Of course, this function already has a series attribute in MuPAD, which you can inspect via expose(exp::series).) In order not to overwrite the already existing attribute, we work on a copy of the exponential function called Exp.

The series attribute must be a procedure with four arguments. This procedure is called whenever a series expansion of Exp with an arbitrary argument is to be computed. The first argument is the argument of Exp in the series call. The second argument is the series variable; the expansion point is always the origin 0; other expansion points are internally moved to the origin by a change of variables. The third and the fourth argument are identical with the order and the dir argument of series, respectively.

For example, the command series(Exp(x^2 + 2), x, 5) is internally converted into the call Exp::series(x^2 + x, x, 5, Undirected). Here is an example of a series attribute for Exp.

// The series attribute for Exp. It handles the call
// series(Exp(f), x = 0, order, dir)
ExpSeries := proc(f, x, order, dir)
  local t, x0, s, r, i;
  // Expand the argument into a series.
  t := series(f, x, order, dir);

  // Determine the order k of the lowest term in t, so that
  // t = c*x^k + higher order terms, for some non-zero c.
  k := ldegree(t);

  if k = FAIL then
    // t consists only of an error term O(..)

  elif k < 0 then
    // This corresponds to an expansion of exp around infinity,
    // which does not exist for the exponential
    // function, since it has an essential singularity. Thus we
    // return FAIL, which makes series return unevaluatedly. For
    // other special functions, you may add an asymptotic
    // expansion here.

  else // k >= 0
    // This corresponds to an expansion of exp around a
    // finite point x0. We write t = x0 + y, where all
    // terms in y have positive order, use the
    // formula exp(x0 + y) = exp(x0)*exp(y) and compute
    // the series expansion of exp(y) as the functional
    // composition of the Taylor series of exp(x) around
    // x = 0 with t - x0. If your special function has
    // any finite singularities, then they should be
    // treated here.
    x0 := coeff(t, x, 0);
    s := Series::Puiseux::create(1, 0, order,
           [1/i! $ i = 0..(order - 1)], x, 0, dir);
    return(Series::Puiseux::scalmult(s @ (t - x0), Exp(x0), 0))

This special function must be embedded in a function environment. The following command defines Exp as a function environment and lets the system function exp do the evaluation. The subs command applied on the result achieves that Exp with symbolic arguments is returned as Exp and not as exp.

Exp := funcenv(x -> subs(exp(x), hold(exp)=hold(Exp))):
Exp(1), Exp(-1.0), Exp(x^2 + x)

series can already handle this “new” function, but it can only compute a Taylor expansion with symbolic derivatives:

ORDER := 3: series(Exp(x), x = 0)

One can define the series attribute of Exp by assigning the procedure above to its series slot:

Exp::series := ExpSeries:

Now we can test the new attribute:

series(Exp(x^2 + x), x = 0) = series(exp(x^2 + x), x = 0)

series(Exp(x^2 + x), x = 2) = series(exp(x^2 + x), x = 2)

series(Exp(x^2 + x), x = 0, 1)

series(Exp(x^2 + x), x = infinity)

Another possibility to obtain series expansions of user-defined functions is to define the diff attribute of the corresponding function environment. This is used by series to compute a Taylor expansion when no series attribute exists. However, this only works when a Taylor expansion exists, whilst a series attribute can handle more general types of series expansions as well.

delete ExpSeries, Exp:



An arithmetical expression representing a function in x


An identifier


The expansion point: an arithmetical expression. If not specified, the default expansion point 0 is used.


The number of terms to be computed: a nonnegative integer or infinity. The default order is given by the environment variable ORDER (default value 6).


Left, Real, Right, Undirected

If no expansion exists that is valid in the complex plane, this argument can be used to request expansions that only need to be valid along the real line. The default is Undirected.


Supresses warning messages printed during the series computation. This can be useful if series is called within user-defined procedures.


Option, specified as UseGseries = b

Use Series::gseries to compute the series. b must be TRUE or FALSE. Default is TRUE. Even if this option is set to TRUE, computing a Puiseux expansion will be attempted first.

Return Values

If order is a nonnegative integer, then series returns either an object of the domain type Series::Puiseux or Series::gseries, an expression of type "series", or, if f is a RootOf expression, a set of type Type::Set. If order = infinity, then series returns an arithmetical expression.

Overloaded By