Numerical solution of an ordinary differential equation

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::odesolve2(f, t0, Y0, <method>, <RememberLast>, <RelativeError = rtol>, <AbsoluteError = atol>, <Stepsize = h>, <MaxStepsize = hmax>)


numeric::odesolve2( f, t0, Y0, … ) returns a function representing the numerical solution Y(t) of the first order differential equation (dynamical system) , Y(t0) = Y0 with and .

The utility function numeric::ode2vectorfield may be used to produce the input parameters f, t0, Y0 from a set of differential expressions representing the ODE. Cf. Example 1.

The function generated by Y := numeric::odesolve2(f, t0, Y0) is essentially

Y := t -> numeric::odesolve(t_0..t, f, Y_0).

Numerical integration is launched, when Y is called with a real numerical argument. The call Y(t) returns the solution vector in a format corresponding to the type of the initial condition Y0 with which Y was defined: Y(t) either yields a list or a 1-dimensional array.

If t is not a real numerical value, then Y(t) returns a symbolic function call.

See the help page of numeric::odesolve for details on the parameters and the options.

The options Alldata = n and Symbolic accepted by numeric::odesolve have no effect: numeric::odesolve2 ignores these options.


Without RememberLast, the function Y remembers all values it has computed. When calling Y(T), it searches its remember table for the time t0 < T < t closest to t, and integrates from T to t using the previously computed Y(T) as initial value. Here, t0 is the time for which the initial value of the ODE is given. This reduces the costs of a call considerably, if Y must be evaluated many times, for example, when plotting the ODE solution. The best approach is to call Y only with a monotonically increasing (or decreasing) sequence of values t starting from t0. Further, the function must be re-initialized whenever DIGITS is increased. See Example 3.

Environment Interactions

The function returned by numeric::odesolve2 is sensitive to the environment variable DIGITS, which determines the numerical working precision.

Without RememberLast, the function returned by numeric::odesolve2 uses option remember.


Example 1

The numerical solution of the initial value problem , y(0) = 2 is represented by the following function Y = [y]:

f := (t, Y) -> [t*sin(Y[1])]:

Alternatively, the utility function numeric::ode2vectorfield can be used to generate the input parameters in a more intuitive way:

[f, t0, Y0] :=
  [numeric::ode2vectorfield({y'(t) = t*sin(y(t)), y(0) = 2}, [y(t)])]

Y := numeric::odesolve2(f, t0, Y0)

The procedure Y starts the numerical integration when called with a numerical argument:

Y(-2), Y(0), Y(0.1), Y(PI + sqrt(2))

Calling Y with a symbolic argument yields a symbolic call:

Y(t), Y(t + 5), Y(t^2 - 4)

eval(subs(%, t = PI))

The numerical solution can be plotted. Note that Y(t) returns a list, so we plot the list element Y(t)[1]:

plotfunc2d(Y(t)[1], t = -5..5):

delete f, t0, Y0, Y:

Example 2

We consider the differential equation with initial conditions y(0) = 0, . The second order equation is converted to a first order system for :


f := (t, Y) -> [Y[2], Y[1]^2]:
t0 := 0: Y0 := [0, 1]:
Y := numeric::odesolve2(f, t0, Y0): 
Y(1), Y(PI)

delete f, t0, Y0, Y:

Example 3

We consider the system


f := (t, Y) -> [Y[1] + Y[2], Y[1] - Y[2]]:
Y := numeric::odesolve2(f, 0, [1, I]):
DIGITS := 5:

Increasing DIGITS does not lead to a more accurate result because of the remember mechanism:

DIGITS := 15: Y(1)

This is the previous value computed with 5 digits, printed with 15 digits. Indeed, only 5 digits are correct. For getting a result that is accurate to full precision, one has to erase the remember table via Y:=subsop(Y,5=NIL). Alternatively, one may create a new numerical solution with a fresh (empty) remember table:

Y := numeric::odesolve2(f, 0, [1, I]): 

delete f, Y, DIGITS:

Example 4

We demonstrate the effect of the option RememberLast. We consider the ODE


f := (t, Y) -> [-Y[1] + sin(t)]:
Y := numeric::odesolve2(f, 0, [1]):
Z := numeric::odesolve2(f, 0, [1], RememberLast):

After many calls of Y, its remember table has grown large. In each call, searching the remember table for input parameters close to the present time value becomes expensive. Created with RememberLast, the procedure Z does not remember all its previously computed values apart from the last one. Consequently, it becomes faster than Y:

time(for i from 1 to 1000 do Y(i/100) end)*msec,
time(for i from 1 to 1000 do Z(i/100) end)*msec

Apart from the efficiency, the values returned by Y and Z coincide:

Y(10.5), Z(10.5)

delete f, Y, Z, i:



A procedure representing the vector field of the dynamical system


A numerical real value for the initial time


A list or 1-dimensional array of numerical values representing the initial value


One of the Runge-Kutta schemes listed below.


BUTCHER6, CK45, CK54, DOPRI45, DOPRI54, DOPRI56, DOPRI65, DOPRI78, DOPRI87, EULER1, GAUSS, RK4, RKF34, RKF43, RKF45a, RKF45b, RKF54a, RKF54b, RKF78, RKF87, xCK45, xCK54, xDOPRI45, xDOPRI54, xDOPRI56, xDOPRI65, xDOPRI78, xDOPRI87, xRKF34, xRKF43, xRKF45a, xRKF45b, xRKF54a, xRKF54b, xRKF78, xRKF87

Option, specified as GAUSS = s

Name of the Runge-Kutta scheme. For details, see the documentation of numeric::odesolve.


Modifies the internal remember mechanism: the procedure returned by numeric::odesolve2 does not remember the results of all previous calls, but only the result of the last call.

Without this option, the procedure returned by numeric::odesolve2 employs option remember to remember the results of all preceding calls. If the function is called very often (hundreds or thousands of times), the remember table grows large and searching this table for entries close to the current time value may become costly. With RememberLast, the procedure returned by numeric::odesolve2 does not use option remember to remember all previous results but implements a very simple and inexpensive mechanism to remember only the result of the very last call.

This option is highly recommended when the numerical procedure returned by numeric::odesolve2 is to be called often (hundreds of thousands of times) with monotonically increasing or decreasing time values. Cf. Example 4.


Option, specified as RelativeError = rtol

Forces internal numerical Runge-Kutta steps to use step sizes with relative local discretization errors below rtol. This tolerance must be a positive numerical real value not smaller than . The default tolerance is RelativeError = 10^(-DIGITS). See the help page of numeric::odesolve for further details.


Option, specified as AbsoluteError = atol

Forces internal numerical Runge-Kutta steps to use step sizes with absolute local discretization errors below atol. This tolerance must be a non-negative numerical real value. The default tolerance is AbsoluteError = 10^(-10*DIGITS). See the help page of numeric::odesolve for further details.


Option, specified as Stepsize = h

Switches off the internal error control and uses a Runge-Kutta iteration with constant step size h. h must be a positive real value. See the help page of numeric::odesolve for further details.


Option, specified as MaxStepsize = hmax

Restricts adaptive step sizes to values not larger than hmax; hmax must be a positive numerical value. See the help page of numeric::odesolve for further details.

Return Values