function ddex3mod
%DDEX3  Example for DDESD, modified to compute the residual.
% 
%   Modifications by Robert M. Corless, January 2011
%   See the DDEX3 example in odedemos['Delay Differential Equations']
%   by Jacek Kierzenka and Lawrence F. Shampine
close all
t0 = 0.1;
tfinal = 5;
tspan = [t0, tfinal];

teensy = 1.0e-5;

% Use a solution structure for the history, to make residuals easy
ddex3hist = ddesd( @(t,y,z) [1/t; -1.0./t.^2], [], ...
                   [log(teensy), 1.0/teensy], [teensy, 0.1] );
% The call is modified to use the structure history
sol = ddesd(@ddex3de,@ddex3delay,ddex3hist,tspan);

% Sample the residual
npts = 1001;
% Step *just* off the initial point (residual huge there)
t = linspace(t0+teensy,tfinal,npts);
[y,yp] = deval( sol, t );
% Compute the corresponding delays
tdelay = ddex3delay( t, y );
% Function at delayed arguments
z = deval( sol, tdelay );
resid = zeros(2,npts);
% ddex3de was not vectorized, step through to compute residuals
for i=1:npts,
    resid(:,i) = yp(:,i) - ddex3de(t(i), y(:,i), z(:,i) );
end
figure
semilogy(t, abs(resid(1,:)),'k')
xlabel('t')
ylabel('residual')
h = diff( sol.x );
figure
semilogy( sol.x(2:end), h, 'k.')
ylabel('step sizes')
xlabel('t')
% The following code is identical to that of DDEX3

function d = ddex3delay(t,y)
% State dependent delay function for DDEX3.  
  d = exp(1 - y(2,:));

% -----------------------------------------------

function dydt = ddex3de(t,y,Z)
% Differential equations function for DDEX3.    
  dydt = [ y(2); -Z(2)*y(2)^2*exp(1 - y(2))];

