Main Content


Solve delay differential equations (DDEs) with general delays


sol = ddesd(ddefun,delays,history,tspan)
sol = ddesd(ddefun,delays,history,tspan,options)



Function handle that evaluates the right side of the differential equations y′(t) = f(t,y(t),y(d(1),...,y(d(k))). The function must have the form

dydt = ddefun(t,y,Z)

where t corresponds to the current t, y is a column vector that approximates y(t), and Z(:,j) approximates y(d(j)) for delay d(j) given as component j of delays(t,y). The output is a column vector corresponding to f(t,y(t),y(d(1),...,y(d(k))).


Function handle that returns a column vector of delays d(j). The delays can depend on both t and y(t). ddesd imposes the requirement that d(j) ≤ t by using min(d(j),t).

If all the delay functions have the form d(j) = t – τj, you can set the argument delays to a constant vector delays(j) = τj. With delay functions of this form, ddesd is used exactly like dde23.


Specify history in one of three ways:

  • A function of t such that y = history(t) returns the solution y(t) for tt0 as a column vector

  • A constant column vector, if y(t) is constant

  • The solution sol from a previous integration, if this call continues that integration


Interval of integration from t0=tspan(1) to tf=tspan(end) with t0 < tf.


Optional integration argument. A structure you create using the ddeset function. See ddeset for details.


sol = ddesd(ddefun,delays,history,tspan) integrates the system of DDEs


on the interval [t0,tf], where delays d(j) can depend on both t and y(t), and t0 < tf. Inputs ddefun and delays are function handles. See Create Function Handle for more information.

Parameterizing Functions explains how to provide additional parameters to the functions ddefun, delays, and history, if necessary.

ddesd returns the solution as a structure sol. Use the auxiliary function deval and the output sol to evaluate the solution at specific points tint in the interval tspan = [t0,tf].

yint = deval(sol,tint)

The structure sol returned by ddesd has the following fields.


Mesh selected by ddesd


Approximation to y(x) at the mesh points in sol.x.


Approximation to y′(x) at the mesh points in sol.x


Solver name, 'ddesd'

sol = ddesd(ddefun,delays,history,tspan,options) solves as above with default integration properties replaced by values in options, an argument created with ddeset. See ddeset and Solving Delay Differential Equations for more information.

Commonly used options are scalar relative error tolerance 'RelTol' (1e-3 by default) and vector of absolute error tolerances 'AbsTol' (all components are 1e-6 by default).

Use the 'Events' option to specify a function that ddesd calls to find where functions g(t,y(t),y(d(1)),...,y(d(k))) vanish. This function must be of the form

[value,isterminal,direction] = events(t,y,Z)

and contain an event function for each event to be tested. For the kth event function in events:

  • value(k) is the value of the kth event function.

  • isterminal(k) = 1 if you want the integration to terminate at a zero of this event function and 0 otherwise.

  • direction(k) = 0 if you want ddesd to compute all zeros of this event function, +1 if only zeros where the event function increases, and -1 if only zeros where the event function decreases.

If you specify the 'Events' option and events are detected, the output structure sol also includes fields:


Row vector of locations of all events, i.e., times when an event function vanished

Matrix whose columns are the solution values corresponding to times in sol.xe

Vector containing indices that specify which event occurred at the corresponding time in sol.xe


The equation

sol = ddesd(@ddex1de,@ddex1delays,@ddex1hist,[0,5]);

solves a DDE on the interval [0,5] with delays specified by the function ddex1delays and differential equations computed by ddex1de. The history is evaluated for t ≤ 0 by the function ddex1hist. The solution is evaluated at 100 equally spaced points in [0,5]:

tint = linspace(0,5);
yint = deval(sol,tint);

and plotted with


This problem involves constant delays. The delay function has the form

function d = ddex1delays(t,y)
%DDEX1DELAYS  Delays for using with DDEX1DE.
d = [ t - 1
      t - 0.2];

The problem can also be solved with the syntax corresponding to constant delays

delays = [1, 0.2];
sol = ddesd(@ddex1de,delays,@ddex1hist,[0, 5]);

or using dde23:

sol = dde23(@ddex1de,delays,@ddex1hist,[0, 5]);

For more examples of solving delay differential equations see ddex2 and ddex3.


[1] Shampine, L.F., “Solving ODEs and DDEs with Residual Control,” Applied Numerical Mathematics, Vol. 52, 2005, pp. 113-127.