This topic shows you how to solve a system of equations symbolically using Symbolic Math Toolbox™. This toolbox offers both numeric and symbolic equation solvers. For a comparison of numeric and symbolic solvers, see Select Numeric or Symbolic Solver.

Suppose you have the system

$$\begin{array}{c}{x}^{2}{y}^{2}=0\\ x-\frac{y}{2}=\alpha ,\end{array}$$

and you want to solve for *x* and *y*.
First, create the necessary symbolic objects.

syms x y a

There are several ways to address the output of `solve`

.
One way is to use a two-output call.

[solx,soly] = solve(x^2*y^2 == 0, x-y/2 == a)

The call returns the following.

solx = 0 a soly = -2*a 0

Modify the first equation to *x*^{2}*y*^{2} = 1. The new system
has more solutions.

[solx,soly] = solve(x^2*y^2 == 1, x-y/2 == a)

Four distinct solutions are produced.

solx = a/2 - (a^2 - 2)^(1/2)/2 a/2 - (a^2 + 2)^(1/2)/2 a/2 + (a^2 - 2)^(1/2)/2 a/2 + (a^2 + 2)^(1/2)/2 soly = - a - (a^2 - 2)^(1/2) - a - (a^2 + 2)^(1/2) (a^2 - 2)^(1/2) - a (a^2 + 2)^(1/2) - a

Since you did not specify the dependent variables, `solve`

uses `symvar`

to
determine the variables.

This way of assigning output from `solve`

is
quite successful for “small” systems. For instance,
if you have a 10-by-10 system of equations, typing the following is
both awkward and time consuming.

[x1,x2,x3,x4,x5,x6,x7,x8,x9,x10] = solve(...)

To circumvent this difficulty, `solve`

can
return a structure whose fields are the solutions. For example, solve
the system of equations `u^2 - v^2 = a^2`

, ```
u
+ v = 1
```

, `a^2 - 2*a = 3`

.

syms u v a S = solve(u^2 - v^2 == a^2, u + v == 1, a^2 - 2*a == 3)

The solver returns its results enclosed in this structure.

S = struct with fields: a: [2×1 sym] u: [2×1 sym] v: [2×1 sym]

The solutions for `a`

reside in the “`a`

-field”
of `S`

.

S.a

ans = -1 3

Similar comments apply to the solutions for `u`

and `v`

.
The structure `S`

can now be manipulated by the field
and index to access a particular portion of the solution. For example,
to examine the second solution, you can use the following statement
to extract the second component of each field.

s2 = [S.a(2), S.u(2), S.v(2)]

s2 = [ 3, 5, -4]

The following statement creates the solution matrix `M`

whose
rows comprise the distinct solutions of the system.

M = [S.a, S.u, S.v]

M = [ -1, 1, 0] [ 3, 5, -4]

Clear `solx`

and `soly`

for
further use.

clear solx soly

Linear systems of equations can also be solved using matrix division. For example, solve this system.

clear u v x y syms u v x y eqns = [x + 2*y == u, 4*x + 5*y == v]; S = solve(eqns); sol = [S.x; S.y] [A,b] = equationsToMatrix(eqns,x,y); z = A\b

sol = (2*v)/3 - (5*u)/3 (4*u)/3 - v/3 z = (2*v)/3 - (5*u)/3 (4*u)/3 - v/3

Thus,`sol`

and `z`

produce
the same solution, although the results are assigned to different
variables.

`solve`

does not automatically return all
solutions of an equation. To return all solutions along with the parameters
in the solution and the conditions on the solution, set the `ReturnConditions`

option
to `true`

.

Consider the following system of equations:

$$\begin{array}{l}\mathrm{sin}\left(x\right)\text{}+\mathrm{cos}\left(y\right)=\frac{4}{5}\\ \mathrm{sin}\left(x\right)\mathrm{cos}\left(y\right)=\frac{1}{10}\end{array}$$

Visualize the system of equations using `fimplicit`

.
To set the *x*-axis and *y*-axis
values in terms of `pi`

, get the axes handles using `axes`

in `a`

.
Create the symbolic array `S`

of the values `-2*pi`

to `2*pi`

at
intervals of `pi/2`

. To set the ticks to `S`

,
use the `XTick`

and `YTick`

properties
of `a`

. To set the labels for the x-and y-axes, convert `S`

to
character vectors. Use `arrayfun`

to apply `char`

to
every element of `S`

to return `T`

.
Set the `XTickLabel`

and `YTickLabel`

properties
of `a`

to `T`

.

syms x y eqn1 = sin(x)+cos(y) == 4/5; eqn2 = sin(x)*cos(y) == 1/10; a = axes; fimplicit(eqn1,[-2*pi 2*pi],'b'); hold on grid on fimplicit(eqn2,[-2*pi 2*pi],'m'); L = sym(-2*pi:pi/2:2*pi); a.XTick = double(L); a.YTick = double(L); M = arrayfun(@char, L, 'UniformOutput', false); a.XTickLabel = M; a.YTickLabel = M; title('Plot of System of Equations') legend('sin(x)+cos(y) == 4/5','sin(x)*cos(y) == 1/10',... 'Location','best','AutoUpdate','off')

The solutions lie at the intersection of the two plots. This
shows the system has repeated, periodic solutions. To solve this system
of equations for the full solution set, use `solve`

and
set the `ReturnConditions`

option to `true`

.

S = solve(eqn1, eqn2, 'ReturnConditions', true)

S = struct with fields: x: [2×1 sym] y: [2×1 sym] parameters: [1×2 sym] conditions: [2×1 sym]

`solve`

returns a structure `S`

with
the fields `S.x`

for the solution to `x`

, `S.y`

for
the solution to `y`

, `S.parameters`

for
the parameters in the solution, and `S.conditions`

for
the conditions on the solution. Elements of the same index in `S.x`

, `S.y`

,
and `S.conditions`

form a solution. Thus, `S.x(1)`

, `S.y(1)`

,
and `S.conditions(1)`

form one solution to the system
of equations. The parameters in `S.parameters`

can
appear in all solutions.

Index into `S`

to return the solutions, parameters,
and conditions.

S.x S.y S.parameters S.conditions

ans = z1 z1 ans = z z ans = [ z, z1] ans = (in((z - acos(6^(1/2)/10 + 2/5))/(2*pi), 'integer') |... in((z + acos(6^(1/2)/10 + 2/5))/(2*pi), 'integer')) &... (in(-(pi - z1 + asin(6^(1/2)/10 - 2/5))/(2*pi), 'integer') |... in((z1 + asin(6^(1/2)/10 - 2/5))/(2*pi), 'integer')) (in((z1 - asin(6^(1/2)/10 + 2/5))/(2*pi), 'integer') |... in((z1 - pi + asin(6^(1/2)/10 + 2/5))/(2*pi), 'integer')) &... (in((z - acos(2/5 - 6^(1/2)/10))/(2*pi), 'integer') |... in((z + acos(2/5 - 6^(1/2)/10))/(2*pi), 'integer'))

To solve the system of equations under conditions, specify the
conditions in the input to `solve`

.

Solve the system of equations considered above for `x`

and `y`

in
the interval `-2*pi`

to `2*pi`

.
Overlay the solutions on the plot using `scatter`

.

Srange = solve(eqn1, eqn2, -2*pi<x, x<2*pi, -2*pi<y, y<2*pi, 'ReturnConditions', true); scatter(Srange.x, Srange.y,'k')

You can use the solutions, parameters, and conditions returned
by `solve`

to find solutions within an interval
or under additional conditions. This section has the same goal as
the previous section, to solve the system of equations within a search
range, but with a different approach. Instead of placing conditions
directly, it shows how to work with the parameters and conditions
returned by `solve`

.

For the full solution `S`

of the system of
equations, find values of `x`

and `y`

in
the interval `-2*pi`

to `2*pi`

by
solving the solutions `S.x`

and `S.y`

for
the parameters `S.parameters`

within that interval
under the condition `S.conditions`

.

Before solving for `x`

and `y`

in
the interval, assume the conditions in `S.conditions`

using `assume`

so
that the solutions returned satisfy the condition. Assume the conditions
for the first solution.

assume(S.conditions(1))

Find the parameters in `S.x`

and `S.y`

.

paramx = intersect(symvar(S.x), S.parameters) paramy = intersect(symvar(S.y), S.parameters)

paramx = z1 paramy = z

Solve the first solution of `x`

for the parameter `paramx`

.

solparamx(1,:) = solve(S.x(1) > -2*pi, S.x(1) < 2*pi, paramx)

solparamx = [ pi + asin(6^(1/2)/10 - 2/5), asin(6^(1/2)/10 - 2/5) - pi, -asin(6^(1/2)/10 - 2/5), - 2*pi - asin(6^(1/2)/10 - 2/5)]

Similarly, solve the first solution of `y`

for `paramy`

.

solparamy(1,:) = solve(S.y(1) > -2*pi, S.y(1) < 2*pi, paramy)

solparamy = [ acos(6^(1/2)/10 + 2/5), acos(6^(1/2)/10 + 2/5) - 2*pi, -acos(6^(1/2)/10 + 2/5), 2*pi - acos(6^(1/2)/10 + 2/5)]

Clear the assumptions set by `S.conditions(1)`

using `assume`

.
Call `asumptions`

to check that the assumptions
are cleared.

assume(S.parameters,'clear') assumptions

ans = Empty sym: 1-by-0

Assume the conditions for the second solution.

assume(S.conditions(2))

Solve the second solution to `x`

and `y`

for
the parameters `paramx`

and `paramy`

.

solparamx(2,:) = solve(S.x(2) > -2*pi, S.x(2) < 2*pi, paramx) solparamy(2,:) = solve(S.y(2) > -2*pi, S.y(2) < 2*pi, paramy)

solparamx = [ pi + asin(6^(1/2)/10 - 2/5), asin(6^(1/2)/10 - 2/5) - pi, -asin(6^(1/2)/10 - 2/5), - 2*pi - asin(6^(1/2)/10 - 2/5)] [ asin(6^(1/2)/10 + 2/5), pi - asin(6^(1/2)/10 + 2/5), asin(6^(1/2)/10 + 2/5) - 2*pi, - pi - asin(6^(1/2)/10 + 2/5)] solparamy = [ acos(6^(1/2)/10 + 2/5), acos(6^(1/2)/10 + 2/5) - 2*pi, -acos(6^(1/2)/10 + 2/5), 2*pi - acos(6^(1/2)/10 + 2/5)] [ acos(2/5 - 6^(1/2)/10), acos(2/5 - 6^(1/2)/10) - 2*pi, -acos(2/5 - 6^(1/2)/10), 2*pi - acos(2/5 - 6^(1/2)/10)]

The first rows of `paramx`

and `paramy`

form
the first solution to the system of equations, and the second rows
form the second solution.

To find the values of `x`

and `y`

for
these values of `paramx`

and `paramy`

,
use `subs`

to substitute for `paramx`

and `paramy`

in `S.x`

and `S.y`

.

solx(1,:) = subs(S.x(1), paramx, solparamx(1,:)); solx(2,:) = subs(S.x(2), paramx, solparamx(2,:)) soly(1,:) = subs(S.y(1), paramy, solparamy(1,:)); soly(2,:) = subs(S.y(2), paramy, solparamy(2,:))

solx = [ pi + asin(6^(1/2)/10 - 2/5), asin(6^(1/2)/10 - 2/5) - pi, -asin(6^(1/2)/10 - 2/5), - 2*pi - asin(6^(1/2)/10 - 2/5)] [ asin(6^(1/2)/10 + 2/5), pi - asin(6^(1/2)/10 + 2/5), asin(6^(1/2)/10 + 2/5) - 2*pi, - pi - asin(6^(1/2)/10 + 2/5)] soly = [ acos(6^(1/2)/10 + 2/5), acos(6^(1/2)/10 + 2/5) - 2*pi, -acos(6^(1/2)/10 + 2/5), 2*pi - acos(6^(1/2)/10 + 2/5)] [ acos(2/5 - 6^(1/2)/10), acos(2/5 - 6^(1/2)/10) - 2*pi, -acos(2/5 - 6^(1/2)/10), 2*pi - acos(2/5 - 6^(1/2)/10)]

Note that `solx`

and `soly`

are
the two sets of solutions to `x`

and to `y`

.
The full sets of solutions to the system of equations are the two
sets of points formed by all possible combinations of the values in `solx`

and `soly`

.

Plot these two sets of points using `scatter`

.
Overlay them on the plot of the equations. As expected, the solutions
appear at the intersection of the plots of the two equations.

for i = 1:length(solx(1,:)) for j = 1:length(soly(1,:)) scatter(solx(1,i), soly(1,j), 'k') scatter(solx(2,i), soly(2,j), 'k') end end

Symbolic calculations provide exact accuracy, while numeric
calculations are approximations. Despite this loss of accuracy, you
might need to convert symbolic results to numeric approximations for
use in numeric calculations. For a high-accuracy conversion, use variable-precision
arithmetic provided by the `vpa`

function.
For standard accuracy and better performance, convert to double precision
using `double`

.

Use `vpa`

to convert the symbolic solutions `solx`

and `soly`

to
numeric form.

vpa(solx) vpa(soly)

ans = [ 2.9859135500977407388300518406219,... -3.2972717570818457380952349259371,... 0.15567910349205249963259154265761,... -6.1275062036875339772926952239014] ... [ 0.70095651347102524787213653614929,... 2.4406361401187679905905068471302,... -5.5822287937085612290531502304097,... -3.8425491670608184863347799194288] ans = [ 0.86983981332387137135918515549046,... -5.4133454938557151055661016110685,... -0.86983981332387137135918515549046,... 5.4133454938557151055661016110685] ... [ 1.4151172233028441195987301489821,... -4.8680680838767423573265566175769,... -1.4151172233028441195987301489821,... 4.8680680838767423573265566175769]

If results look complicated, `solve`

is stuck,
or if you want to improve performance, see, Troubleshoot Equation Solutions from solve Function.