function x = dde(maxDelay,initialConditions,timeRange,ddeFunction)
    // a solver similar to ode for delayed differential equations
    // 11/03/2012 : Philippe Thieullen
    // maxDelay = scalar : 
    //     the maximum delay, each specific delay tau_i used in 
    //     ddeFunction should be in the range 0 < tau_i <= maxDelay
    // initialConditions = function : 
    //     a function that returns the initial backward conditions 
    //     for different bt(i) <= t_init of the form
    //     function  y = initialCondition(bt)
    //         bt = row vector of size r : backward times
    //         y = matrix of size n*r : the backward positions
    //         y(:,i) = initialConditions(bt(i))
    //     endfunction
    // timeRange = row vector of size T : 
    //     a set of observation times [t_init,...,t_final] 
    //     with t_init < t_final
    // ddeFunction = function : 
    //     the differential equation of the form
    //     function xdot = ddeFunction(t,x)
    //         t = scalar :the present time
    //         x = column vector of size n : positions x(t) at time t
    //         xdot = column vector of size n : the derivative dx/dt
    //     endfunction
    //     that solves the equation 
    //     dx/dt = y = f(t,x(t),x(t-tau_1),...,x(t-tau_d))
    //     the backward positions z(:,i) = x(t-tau_i) are obtained by
    //     z = ddeGetBackwardPositions(t,x,t-[tau_1,...,tau_d])
    // x = matrix of size n*T : 
    //     the positions x(:,i) at timeRange(i)
    
    global ddeDelayOrigin;
    global ddeDelayTimes;
    global ddeBackwardPositions;

    delaysNumber     = 1000;
    ddeDelayOrigin   = timeRange(1);
    ddeSecurityDelay = 1;
    ddeDelayTimes = linspace(-maxDelay-ddeSecurityDelay,0,delaysNumber);
    ddeBackwardPositions = initialConditions(ddeDelayOrigin + ddeDelayTimes);
    x = ode('stiff',ddeBackwardPositions($), ddeDelayOrigin,timeRange,ddeFunction);
endfunction

function z = ddeGetBackwardPositions(t,x,bt)
    // a private function to dde
    // t = scalar : 
    //     the present time
    // x = column vector of size n : 
    //     the positions at time t
    // bt = row vector of size d : 
    //     backward times that should satisty t-maxDelay <= bt(i) <= t
    // z = matrix of size n*d : 
    //     z(:,i) = x(bt(i)) the positions at time bt(i)

    global ddeDelayOrigin;
    global ddeDelayTimes;
    global ddeBackwardPositions;
    
    timesScale = ddeDelayOrigin + ddeDelayTimes;
    positionsScale = ddeBackwardPositions;
    // t is not necessarily increasing
    if t > ddeDelayOrigin then 
        timesScale = [timesScale,t];
        positionsScale = [ddeBackwardPositions,x];
        ddeBackwardPositions = ...
            interp1(timesScale, positionsScale, t + ddeDelayTimes);
        ddeDelayOrigin = t; 
    end 
    z = interp1(timesScale, positionsScale, bt); 
endfunction
