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
%-------------------------------------------
