System of state dependent delay equations causes ddesd to run forever.

I am trying to get a numerical solution to a threshold-type state dependent delay equation I have been working on. However, it causes the solver ddesd to loop forever. I am unfamiliar with the numerical method that ddesd uses, so any help would be appreciated.

Here are the equations for the system:

function dydt = ddefun(t,y,Z) % equation being solved
kplus = 5;
kminus = 0.1;
a = 5;
l = 25;
dydt = [a*y(4) - (kplus*y(1) - kminus)*(y(2) + y(3));
    (kplus*y(1) - kminus)*(y(1) - Z(1,1));
    (kplus*y(1) - kminus)*Z(1,1) - a*y(3);
    (kplus*y(1) - kminus)*(l*Z(1,1) + y(3)) - a*y(4);
    1-((kplus*y(1) - kminus)/(kplus*Z(1,1) - kminus))];
end
%-------------------------------------------
function d = dely(t,y) % delay for y
d = t - y(5);
end
%-------------------------------------------
function v = history(t) % history function for t < t0
v = [20;
    15;
    0;
    0;
    0;];
end
%-------------------------------------------

When I call

tspan = [0 20];
sol = ddesd(@ddefun, @dely, @history, tspan);

the script runs without error, but doesn't produce a result. I am aware that numerical methods can take some time, but I have left the script running for over 48 hours with no solution. When I run ddesd on an extremely similar system (for example, remove the a*y(4) term in the first equation) the script finishes almost instantly. This is why I believe the solver is stuck in a loop instead of simply taking a long time.

Answer 1

[Torsten]

y(5) is negative shortly after t = 0. Thus your delay d = t - y(5) in principle would refer to values of y(1) that are not yet known.

Ok, "ddesd" uses min(t,d) as delay, thus Z(1,1) = y1(t), but neverthess for me it seems that your problem is not well-posed.

tspan = [0 1.5];
options = ddeset('RelTol',1e-8,'AbsTol',1e-8);
sol = ddesd(@ddefun, @dely, @history, tspan, options);
plot(sol.x,sol.y(5,:))

function dydt = ddefun(t,y,Z) % equation being solved
kplus = 5;
kminus = 0.1;
a = 5;
l = 25;
dydt = [a*y(4) - (kplus*y(1) - kminus)*(y(2) + y(3));
    (kplus*y(1) - kminus)*(y(1) - Z(1,1));
    (kplus*y(1) - kminus)*Z(1,1) - a*y(3);
    (kplus*y(1) - kminus)*(l*Z(1,1) + y(3)) - a*y(4);
    1-((kplus*y(1) - kminus)/(kplus*Z(1,1) - kminus))];
end
%-------------------------------------------
function d = dely(t,y) % delay for y
d = t - y(5);
end
%-------------------------------------------
function v = history(t) % history function for t < t0
v = [20;
    15;
    0;
    0;
    0;];
end
%-------------------------------------------