Solve Ordinary Differential Equations and Systems

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.

General Solutions

An ordinary differential equation (ODE) contains derivatives of dependent variables with respect to the only independent variable. If y is a dependent variable and x is an independent variable, the solution of an ODE is an expression y(x). The order of the derivative of a dependent variable defines the order of an ODE.

The solution of a single explicit first-order ODE can always be computed by integration, provided the solution exists. To define an ordinary differential equation, use the ode command:

o := ode(y'(x) = y(x)^2, y(x))


ode does not accept multivariate expressions such as y(x, t). It also does not accept piecewise expressions.

Now use the general solve to solve this equation:


Alternatively, you can call the ODE solver directly:

ode::solve(y'(x) = y(x)^2, y(x))

The general solutions of ODEs contain arbitrary constants of integration. The solver generates the constants of integration using the format of an uppercase letter C followed by an automatically generated number. For example, it generates C1, C2, and so on.

For higher-order ODEs, you can find explicit solutions only for special types of equations. For example, the following second-order equation has a solution in terms of elementary functions:

ode::solve(y''(x) = y(x), y(x))

The solver introduces the solution of the following equation in terms of the Bessel functions. MuPAD® uses standard mathematical notations for Bessel and other special functions:

ode::solve(y''(x) = y'(x) + y(x)*exp(x), y(x))

If you have a second- or higher-order ODE, nonlinear ODE, or a system of ODEs, a symbolic solution does not always exist:

ode::solve(y''(x) = y'(x)^2 + y(x)*exp(x), y(x))

For ODEs that cannot be solved symbolically, try using numeric solvers.

Initial and Boundary Value Problems

Many problems in engineering and physics involve solving differential equations with initial conditions or boundary conditions or both. To specify initial or boundary conditions, create a set containing an equation and conditions. For example, state the following initial value problem by defining an ODE with initial conditions:

IVP := ode({y''(x) = y(x), y(0) = 5, y'(0) = 1}, y(x))

When you solve an ODE with initial or boundary conditions, the solver adjusts the integration constants to fit these conditions:


The following equation has both initial and boundary conditions:

ode::solve({y'''(x) = y(x), y(0) = 0, y(5) = 1, y'(0) = 0}, y(x))

Each independent condition removes one integration constant:

ode::solve({y'''(x) = y'(x)}, y(x))

ode::solve({y'''(x) = y'(x), y(0) = 0}, y(x))

ode::solve({y'''(x) = y'(x), y(0) = 0, y(1) = 1}, y(x))

ode::solve({y'''(x)=y'(x), y(0)=0, y(1)=1, y'(0)=1/2}, y(x))

Special Types of Ordinary Differential Equations

Suppose, the equation you want to solve belongs to the Clairaut type:

o:= ode(y(x) = x*y'(x) + y'(x)^3, y(x)):

The solver recognizes the type of the equation and applies the algorithm for solving Clairaut equations. To improve performance, call the solver with the option Type = Clairaut:

solve(o, Type = Clairaut)

The solver tries to recognize and tries to solve the following classes of ODEs.

TypeEquationODE Solver Option
Abel differential equationAbel
Bernoulli differential equation where n ≠ 0 and n ≠ 1Bernoulli
Chini differential equationChini
Clairaut differential equationClairaut
Exact first-order ordinary differential equation that can be represented as M(x, y) dx + N(x, y) dy = 0 where ExactFirstOrder
Exact second-order ordinary differential equation, if a first integral turns this equation into a first-order ODEExactSecondOrder
Linear homogeneous ordinary differential equationDy = 0, where D is a linear differential operatorHomogeneous
Lagrange differential equationLagrange
Riccati differential equationRiccati

If the solver cannot identify the equation with the type you indicated, it issues a warning and returns the special value FAIL:

ode::solve(y'(x) + y(x) = x, y(x), Type = Homogeneous)
Warning: Unable to detect homogeneous ODE. [ode::homogeneous]

Systems of Ordinary Differential Equations

To solve a system of differential equations, specify the system as a set of equations:

s := {y'(x) = z(x), z'(x) = y(x) + 2*z(x)}:

Call the ode::solve function specifying the set of functions {y(x), z(x)} for which you want to solve the system:

ode::solve(s, {y(x), z(x)})

Now, suppose the system of differential equations appears in a matrix form. For example, define the system , where A, B, and Y represent the following matrices:

Y:= matrix([x(t), y(t)]):
A:= matrix([[1, 2], [-1, 1]]):
B:= matrix([1, t]):

The ode::solve function does not accept matrices. To be able to use this solver, extract the components of the matrix and include them in a set. Use the op function to extract the equations from the matrix. Then, use the braces to create a set of the equations. You can omit the right sides of equations, in which case MuPAD assumes them to be 0:

s := {op(diff(Y, t) - A*Y - B)}

Now, specify the set of functions {x(t), y(t)} for which you want to solve the system. Solve the system:

ode::solve(s, {x(t), y(t)})

If you are solving several similar systems of ordinary differential equations in a matrix form, create your own solver for these systems, and then use it as a shortcut. The solver for such systems must be a function that accepts matrices as input arguments, and then performs all required steps. For example, create a solver for a system of the first-order linear differential equations in a matrix form , where the components of functions depend on the variable t:

solveLinearSystem := (A, B, Y) ->
solve(ode({op(diff(Y, t) - A*Y - B)}, {op(Y)})):

The solveLinearSystem function accepts matrices as input parameters, creates a matrix of equations, extracts these equations to a set, and solves the system:

Y:= matrix([x(t), y(t)]):
A:= matrix([[1, 2], [-3, 1]]):
B:= matrix([2, t]):
solveLinearSystem(A, B, Y)

Plot Solutions of Differential Equations

Suppose you want to solve the following equation. The solver returns the results as a set even if the set contains only one element:

f := ode::solve({y''(x) = x*y'(x), y(0) = 0, y'(0) = 1/3}, y(x))

The plotting functions in MuPAD do not accept sets. To plot the solution, access the elements of a solution set using square brackets or the op command:

plotfunc2d(f[1], x = -3..3)

If you have more than one element of a solution set, you can access a particular element. For example, pick the second element of the solution set for the system of ODEs:

f := ode::solve(
               {y'(x) = z(x), z'(x) = y(x) + 2*z(x),
                y(0) = 0, z(0) = 1}, {y(x), z(x)})

The solver returns results for a system as a set that contains a list. To open the set and access the list, use square brackets or the op command. To access a particular entry of this list, use square brackets:


To access the right side of the equation, use square brackets or the rhs command:


Plot this solution:

plotfunc2d(f[1][2][2], x = -1..1)

To plot a solution of the system of ODEs in 3-D, use the plot::Curve3d command:

solution := plot::Curve3d([x, f[1][2][2], f[1][1][2]],
                              x = -2..2, GridVisible):

MuPAD provides the functions plot::Ode2d and plot::Ode3d for visualizing solutions of ODEs. Also, you can plot a vector field associated with an ODE. For all graphic capabilities available in MuPAD, see Graphics and Animations.