This example shows how to use ddensd to solve the neutral DDE presented by Paul  for 0 ≤ t ≤ π.
The equation is
y '(t) = 1 + y(t) – 2y(t/2)2 – y '(t – π)
y(t) = cos (t) for t ≤ 0.
Create a new program file in the editor. This file will contain a main function and four local functions.
Define the first-order DDE as a local function.
function yp = ddefun(t,y,ydel,ypdel) yp = 1 + y - 2*ydel^2 - ypdel; end
Define the solution delay as a local function.
function dy = dely(t,y) dy = t/2; end
Define the derivative delay as a local function.
function dyp = delyp(t,y) dyp = t-pi; end
Define the solution history as a local function.
function y = history(t) y = cos(t); end
Define the interval of integration and solve the DDE using the ddensd function. Add this code to the main function.
tspan = [0 pi]; sol = ddensd(@ddefun,@dely,@delyp,@history,tspan);
Evaluate the solution at 100 equally spaced points between 0 and π. Add this code to the main function.
tn = linspace(0,pi); yn = deval(sol,tn);
Plot the results. Add this code to the main function.
figure plot(tn,yn); xlim([0 pi]); ylim([-1.2 1.2]) xlabel('time t'); ylabel('solution y'); title('Example of Paul with 1 equation and 2 delay functions')
Run your program to calculate the solution and display the plot.
 Paul, C.A.H. "A Test Set of Functional Differential Equations." Numerical Analysis Reports. No. 243. Manchester, UK: Math Department, University of Manchester, 1994.