octave and delay ode's

[franco basaglia]

	Hi all,

I'm trying to solve a system of Seven ODE's (from Stella) by Octave using lsode. 
Two of this are delay differential equations. For example on state variable S(t):

%S(t)
% S(t)=S(t−dt)+(dSi−dSo)×dt

% dSi = 0.2×DELAY((F+0.5×U)×(1−GP),5)×(1−S)×(1−T)×AE %this is Stella definition
% DELAY on S(t)
%%% F=x(2); U=x(4); GP=1-exp(d*x(5)); S=x(6); T=x(7)

if (t < 5*delta_t)
t_back = t;
t = 0;
delay_si = ((x(2)+0.5*x(4))*(exp(d*x(5))));
t = t_back;
t;
else
t_back = t;
t = t-5*delta_t;
delay_si = ((x(2)+0.5*x(4))*(exp(d*x(5))));
t = t_back;
endif

dSi = 0.2*delay_si*(1-x(6))*(1-x(7))*AE;

dSo = (0.1 + 1 - exp(d*x(5))) * (0.1 + x(7)) * x(6) * 0.2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xdot(6) = dSi - dSo ;


Using lsode(adaptive step size method), for Cycle doesn't work because there is an “incompatibility” on delta_t.
Maybe if I print a matrix with x(t) value for every dt I can run for Cycle on these values. 
But I don't know how to print this matrix with lsode.

Any idea ?
Alternative way to solve ode's delay?

Thanks

best regards
f.b.