Search for a numerical root of a system of equations

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.


numeric::fsolve(eq, x, options)
numeric::fsolve(eq, x = a, options)
numeric::fsolve(eq, x = a .. b, options)
numeric::fsolve(eqs, [x1, x2, …], options)
numeric::fsolve(eqs, {x1, x2, …}, options)
numeric::fsolve(eqs, [x1 = a1, x2 = a2, …], options)
numeric::fsolve(eqs, {x1 = a1, x2 = a2, …}, options)
numeric::fsolve(eqs, [x1 = a1 .. b1, x2 = a2 .. b2, …], options)
numeric::fsolve(eqs, {x1 = a1 .. b1, x2 = a2 .. b2, …}, options)


numeric::fsolve(eqs, ...) returns a numerical approximation of a solution of the system of equations eqs.

This is the MuPAD® numerical solver for non-linear systems of equations.


By default, this routine returns only one numerical solution!

The equations must not contain symbolic objects other than the unknowns that cannot be converted to numerical values via float. Symbolic objects such as π or etc. are accepted. The same holds true for starting values and search ranges. Search ranges may contain . Cf. Example 2.

numeric::fsolve implements a purely numerical Newton type root search with a working precision set by the environment variable DIGITS. Well separated simple roots should be exact within this precision. However, multiple roots or badly separated roots may be computed with a restricted precision. Cf. Example 3.


For systems of equations, the expressions defining the equations must have a symbolic derivative!

Overdetermined systems (i.e., more equations than indeterminates) are not accepted. However, there may be more indeterminates than equations. Cf. Example 4.

Specifying indeterminates [x1, x2, …] without starting values or search ranges is equivalent to the search ranges [x_1 = -infinity .. infinity, x_2 = -infinity .. infinity, dots]. Note, however, that the user should assist numeric::fsolve by providing specific search ranges whenever possible! If a complex starting point or a search range involving a complex number is specified for at least one of the unknowns, the search is extended to the entire complex plane for all variables for which no explicit search interval is given.

For real equations and real starting points or search ranges, the internal Newton iteration will usually produce real values, i.e., numeric::fsolve searches for real roots only (unless square roots, logarithms etc. happen to produce complex values from real input). Use complex starting points or search ranges to search for complex roots of real equations. Cf. Example 5.

Starting values and search ranges can be mixed. Cf. Example 6.

Search ranges should only be provided if a solution is known to exist inside the search range. Otherwise, the search may take some time before numeric::fsolve gives up.

Specification of a search range primarily means that starting points from this range are used for the internal Newton search. For sufficiently small search ranges enclosing a solution the search will usually pick out this solution. However, it may also happen that the Newton iteration drifts towards other solutions.

With the default search strategy RestrictedSearch, only solutions from the search range are accepted, even if solutions outside the search range are found internally. More specifically, if a search range such as x = a .. b is specified for the variable x, then solutions satisfying min(ℜ(a), ℜ(b)) ≤ ℜ(x) ≤ max(ℜ(a), ℜ(b)) and min(ℑ(a), ℑ(b)) ≤ ℑ(x) ≤ max(ℑ(a), ℑ(b)) are searched for. Thus, the values a, b specify the bottom left and top right corner of a rectangular search area in the complex plane when the RestrictedSearch strategy is used.

With the search strategy UnrestrictedSearch, any solution inside or outside the search range is accepted and returned. Cf. Example 7.

If starting values for all indeterminates are provided, then a single Newton iteration with these initial data is launched. It either leads to a solution or numeric::fsolve gives up and returns FAIL. The same holds true if search ranges x = a .. a or [x_1 = a_1 .. a_1, x_2 = a_2 .. a_2, dots] of zero length are specified.


The risk of failure is high when providing bad starting values! Starting values are appropriate only if a sufficiently good approximation of the solution is known! On the other hand, providing good starting values is the fastest way to a solution. Cf. Example 8.

If at least one of the indeterminates has a non-trivial search range, then numeric::fsolve uses several Newton iterations with different starting values from the search range. Cf. Example 9. Search ranges in conjunction with the option UnrestrictedSearch provide a higher chance of detecting roots than (bad) starting values!


User defined assumptions such as assume(x > 0) are not taken into account in the numerical search! Provide search ranges instead! Cf. Example 2.


Convergence may be slow for multiple roots! Furthermore, numeric::fsolve may fail to detect such roots!

Use linsolve or numeric::linsolve for systems of linear equations.

Use numeric::realroots, if all real roots of a single non-polynomial real equation in a finite range are desired.

Use polylib::realroots, if all real roots of a real univariate polynomial are desired.

Use numeric::polyroots, if all real and complex roots of a univariate polynomial are desired.

Use numeric::solve, if all roots of a multivariate polynomial system are desired.

The routine numeric::solve provides a common interface to all these numerical solvers.

Environment Interactions

The function is sensitive to the environment variable DIGITS, which determines the numerical working precision.


Example 1

We compute roots of the sine function:

numeric::fsolve(sin(x) = 0, x)

With the option Random, several calls may result in different roots:

numeric::fsolve(sin(x), x, Random)

numeric::fsolve(sin(x), x, Random)

Particular solutions can be chosen by an appropriate starting point close to the wanted solution, or by a search interval:

numeric::fsolve(sin(x), x = 3), 
numeric::fsolve(sin(x), x = -4 .. -3)

The solutions found by numeric::fsolve can be used in subs and assign to substitute or assign the indeterminates:

eqs := [x^2 = sin(y), y^2 = cos(x)]:
solution := numeric::fsolve(eqs, [x, y])

eval(subs(eqs, solution))

assign(solution): x, y

delete eqs, solution, x, y:

Example 2

We demonstrate the use of search ranges. The following system has solutions with positive and negative x. The solution with x ≥ 0 is obtained with the search interval x = 0 .. infinity:

numeric::fsolve([x^2 = exp(x*y), x^2 = y^2], [x = 0 .. infinity, y])

We search for a solution with x ≤ 0:

numeric::fsolve([x^2 = exp(x*y), x^2 = y^2], [x = -infinity .. 0, y])

Example 3

Multiple roots can only be computed with a restricted precision:

numeric::fsolve(expand((x - 1/3)^5), x = 0.3)

Example 4

The following system of equations is degenerate and has a 1-parameter family of solutions. Each call to numeric::fsolve picks out one random solution:

numeric::fsolve([x^2 - y^2, x^2 - y^2], [x, y], Random) $ i = 1 .. 3

The equation may also be specified as an underdetermined system:

numeric::fsolve([x^2 - y^2], [x, y])

Example 5

The following equation has no real solution. Consequently, the numerical search with real starting values fails:

numeric::fsolve(sin(x) + cos(x)^2 = 3, x)

With a complex starting value, a solution is found:

numeric::fsolve(sin(x) + cos(x)^2 = 3, x = I)

Also complex search ranges may be specified. In the following, the internal starting point is a random value on the line from 2 + I to 3 + 2*I. Solutions are accepted if they lie in the complex rectangle with the bottom left corner 2 + I and the top right corner 3 + 2*I:

numeric::fsolve(sin(x) + cos(x)^2 = 3, x = 2 + I .. 3 + 2*I)

Example 6

Starting values and search intervals can be mixed:

numeric::fsolve([x^2 + y^2 = 1, y^2 + z^2 = 1, x^2 + z^2 = 1],
                [x = 1, y = 0 .. 10, z])

Example 7

With UnrestrictedSearch, search intervals are only used for choosing starting values for the internal Newton search. The numerical iteration may drift towards a solution outside the search range:

eqs := [x*sin(10*x) = y^3, y^2 = exp(-2*x/3)]: 
numeric::fsolve(eqs, [x = 0 .. 1, y = -1 .. 0], UnrestrictedSearch)

With the default strategy RestrictedSearch, only solutions inside the search range are accepted:

numeric::fsolve(eqs, [x = 0 .. 1, y = -1 .. 0])

In the last search, also the previous solution outside the search range was found. With the option MultiSolutions, numeric::fsolve returns a sequence of all solutions that were found in the internal search:

numeric::fsolve(eqs, [x = 0 .. 1, y = -1 .. 0], MultiSolutions)

delete eqs:

Example 8

Usually, most of the time is spent internally searching for some (crude) approximations of the root. If high precision roots are required, it is recommended to compute first approximations with moderate values of DIGITS and use them as starting values for a refined search:

eq := exp(-x) = x:
DIGITS := 10: 
firstApprox := numeric::fsolve(eq, x)

This output is suitable as input defining a starting value for x:

DIGITS := 1000: numeric::fsolve(eq, firstApprox)

delete eq, firstApprox, DIGITS:

Example 9

Specifying starting values for the indeterminates launches a single Newton iteration. This may fail, if the starting values are not sufficiently close to the solution:

eq := [x*y = x + y - 4, x/y = x - y + 4]:
numeric::fsolve(eq, [x = 1, y = 1])

If a search range is specified for at least one of the unknowns, then several Newton iterations with random starting values in the search range are used, until a solution is found or until numeric::fsolve gives up:

numeric::fsolve(eq, [x = 1, y = 0 .. 10])

delete eq:



An arithmetical expression or an equation in one indeterminate x. An expression eq is interpreted as the equation eq = 0.


A list, set, array, or matrix (Cat::Matrix) of expressions or equations in several indeterminates x1, x2, ... Expressions are interpreted as homogeneous equations.

x, x1, x2, …

Identifiers or indexed identifiers to be solved for.

a, a1, a2, …

Real or complex numerical starting values for the internal search. Typically, crude approximations of solution.

a .. b, a1 .. b1, a2 .. b2, …

Ranges of numerical values defining search intervals for the numerical root.



Makes numeric::fsolve return only numerical roots in the user-defined search range x = a .. b and [x_1 = a_1 .. b_1, x_2 = a_2 .. b_2, Symbol::hellip], respectively. This is the default search strategy, if a search range is specified for at least one of the unknowns.

Once a root with components (r1, r2, …) is found, it is checked whether min((ai), (bi)) ≤ (ri) ≤ max((ai), (bi)) and min((ai), (bi)) ≤ (ri) ≤ max((ai), (bi)) is satisfied. If the root is not inside the search range, the search is continued. Note that solutions outside the search range may be found internally. These may be accessed with the option MultiSolutions. See Example 7.


Allows numeric::fsolve to find and return solutions outside the specified search range. With this option, the search range is only used to choose random starting points for the internal numerical search.

This option switches off the search strategy RestrictedSearch. With UnrestrictedSearch, numeric::fsolve stops its internal search whenever a root is found, even if the root is not inside the specified search range. Starting points for the internal Newton search are taken from the search range.


Makes numeric::fsolve return all solutions found in the internal search

This option only has an effect when used with the default search strategy RestrictedSearch. A sequence of all roots found in the internal search is returned. Cf. Example 7.


With this option, several calls to numeric::fsolve with the same input parameters may produce different roots.

With this option, random starting values are chosen for the internal search. Consequently, calling numeric::fsolve several times with the same parameters may lead to different solutions. This may be useful when several roots of one and the same equation or set of equations are desired.

Return Values

Single numerical root is returned as a list of equations [x = value] or [x1 = value1, x2 = value2, …], respectively. FAIL is returned if no solution is found. With the option MultiSolutions, sequences of solutions may be returned.


Internally the set of equations f(x) = 0 is solved by a modified Newton iteration with some adaptively chosen step size t. For degenerate or ill-conditioned Jacobians a minimization strategy for is implemented. For scalar real equations, numeric::realroot is used, if a real finite search range is specified.