# PowerSeriesODESolver
# ====================
#
# This package extends a UPS Domain with the procedure ode to solve
# regular linear or non-linear ordinary differential equations
# of arbitrary order.


# Interface:

   # PROCEDURE ode (P: UPSDomain; 
   #                f: LIST OF UPS &-> UPS; 
   #                C: LIST OF UPSDomain[CoefficientRing]): UPS;
   #
   # ode(P, f, C) is the solution to
   #          (n)                    (n-1)
   #         y   =f(y,y',y'', ..., y      )
   # such that
   #          (i)
   #         y   =C[i], for i in 0..n-1
   #

# Implementation

macro( odeprod = `ODESolve/ODEprod` );
ODESolve := proc(P, f, C) local lorder, R;

   if not hasCategory(P, UnivariatePowerSeries) then
      ERROR(`first argument must be a UnivariatePowerSeries category`)
   fi;
   R := P[CoefficientRing];
   if not hasCategory(R, Field) then
      ERROR(`coefficient ring must be a field`)
   fi;
   if not hasOperation(P,FixedPoint) then
      ERROR(`this implementation requires a FixedPoint operation`)
   fi;

   for lorder to nops(C) while R[`=`](C[lorder],R[0]) do od;
   P[FixedPoint](
      subs(['_P'=P, '_f'=f, '_C'=C],
         y -> odeprod(_P, _f, _C, y, nops(_C))),
      0,
      infinity,
      lorder-1
   )
end:

# defOperation(odeprod,  [[UPS,Nary(UPS)] &-> UPS,
#            List(R), UPS, Integer] &-> UPS, UPS);
# PROCEDURE odeprod(P: UPSDomain;
#                   f: LIST OF UPS &-> UPS;
#                   C: LIST OF UPSDomain[CoefficientRing];
#                   y: UPS;
#                   n: INTEGER): UPS;
#
odeprod := proc(P, f, C, y, n) local argseq, yp, i;
   if C=[] then
      argseq := y; yp := y;
      for i to n-1 do
         yp := P[Diff](yp);
         argseq := argseq, yp;
      od;
      f(argseq);
   else
      P[`+`](
         P[Constant](C[1]),
         P[Integrate](
            odeprod(P, f, [C[2..nops(C)]], y, n)
         )
      )
   fi;
end:

save ODESolve,
     odeprod,
     
     `PSODES.m`;

quit
