#
# DDE Solver version 0.1
# RMC August 2004
#
# Calling Sequence
# (sol,dsol) := March( N, f, History )
# f = (t,y,z,j) -> dde, eg (theta,y,z,j) -> D(y)(theta) - a*z(theta) - b*y(theta)
#                      means y'(theta) = a*y(theta-1) + b*y(theta)
#                      to refer to t in coeffs use t = j+theta, 0 <= theta <= 1
#
# Thx Edgardo for "useint"
#
March := proc( N::posint, f, History)
   local j, dy, y, Y, theta, sol, dsol;
   y     := array(0..N); # solution components (operators)
   dy    := array(0..N); # derivative components (operators)
   y[0]  := History;
   dy[0] := D(History);
   for j to N do
      y[j] := unapply(
                subs(
                  dsolve( {f(theta,Y,y[j-1],j), Y(0)=y[j-1](1)}, Y(theta), useint ),
                  Y(theta)
                     ),
                theta
                      );
      dy[j] := D(y[j]);
   end do;
   sol  := t ->  y[floor(t+1)](frac(t+1));
   dsol := t -> dy[floor(t+1)](frac(t+1));
   return sol, dsol;
end proc;