reduce-historical/r38/packages/odesolve/odelin.red

module odelin$  % Simple linear ODE solver

% F.J.Wright@Maths.QMW.ac.uk, Time-stamp: <14 September 2000>

% Based in part on code by Malcolm MacCallum.

% TO DO:

%   Polynomial solutions for polynomial coeffs.

%   Check the distinction between finding a solution technique for an
%   ODE that then fails internally to solve it (e.g. failing to solve
%   an auxiliary equation) and failing to find a solution technique.
%   (Should the former be handled by `odefailure'?)


%% Techniques implemented
%% ======================
%%   First order (integrating factor)
%%   Constant coefficients
%%   Euler and shifted Euler
%%   Exact
%%   Trivial order reduction (dep var and low-order derivs missing)

%%   Second-order special function ODEs (module odespcfn)

%% Notes: Overall factors are handled in most cases by making the ODE
%% "monic".


%% Internal representation
%% =======================

%% A linear ode is represented by its list of coefficient functions
%% (odecoeffs), driver term (driver), and dependent (y) and
%% independent (x) variables.  The maximum (ode_order) and minimum
%% (min_order) derivative orders are included in the representation
%% for efficiency/convenience.  Its solution is represented as a basis
%% for the solution space of the reduced ODE together with a
%% particular integral of the full ODE.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% algebraic procedure odesolve(ode, y, x);
%%    %% Temporary definition for test purposes.
%%    begin scalar !*precise, solution;
%%       ode := num !*eqn2a ode;           % returns ode as expression
%%       if (solution := ODESolve!-linear(ode, y, x)) then
%%          return solution
%%       else
%%          write "***** ODESolve cannot solve this ODE!"
%%    end$

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global '(ODESolve_Before_Lin_Hook ODESolve_After_Lin_Hook)$

algebraic procedure ODESolve!-linear(ode, y, x);
   %% MAIN LINEAR SOLVER
   %% Assumes ODE is an algebraically irreducible "polynomial" expression.
   begin scalar reduced_ode, auxvar, auxeqn, odecoeffs,
         first_arb, solution, driver;
      %% The following decomposition is needed FOR ALL linear ODEs.
      %% The DRIVER is the part of the ODE independent of y, such that
      %% the ODE can be expressed as REDUCED_ODE = DRIVER.
%%       driver := if part(ode, 0) = plus then
%%          -select(~u freeof y, ode) else 0;
      driver := -sub(y=0, ode);
      reduced_ode := ode + driver;
      auxvar := symbolic gensym();
      %% df(y, x, n) => m^n, where m = auxvar
      auxeqn := sub(y=e^(auxvar*x), reduced_ode)/e^(auxvar*x);
      odecoeffs := coeff(auxeqn, auxvar); % low .. high
      traceode "This is a linear ODE of order ", high_pow, ".";
      first_arb := !!arbconst + 1;
      symbolic if not(solution := ODESolve!-linear!-basis
         (odecoeffs, driver, y, x, high_pow, low_pow) or
         (not !*odesolve_fast and
         %% Add a switch to control access to this thread?
         %% It is currently necessary for Zimmer (8).
         << traceode
            "But ODESolve cannot solve it using linear techniques, so ...";
            %% NB: This will probably produce a NONLINEAR ODE!
            %% But, in desperation, try it anyway ...
            (ODESolve!-Interchange(ode, y, x) where !*odesolve_basis = nil)
         >>)) then return;
      %% Return solution as BASIS or LINEAR COMBINATION, assuming a
      %% SINGLE solution since the ODE is linear:
      return if symbolic !*odesolve_basis then % return basis
         if (part(solution, 1, 0) = equal) then
            if lhs first solution = y and % solution is explicit
         (auxeqn := ODESolve!-LinComb2Basis
            (rhs first solution, first_arb, !!arbconst)) then auxeqn
            else << write
               "***** Cannot convert nonlinear combination solution to basis!";
            solution >>
         else solution
      else                              % return linear combination
         if part(solution, 1, 0) = list then
            {y = ODESolve!-Basis2LinComb solution}
         else solution
   end$

algebraic procedure ODESolve!-linear!-basis
   (odecoeffs, driver, y, x, ode_order, min_order);
   %% Always returns the solution in basis format.
   %% Called by ODESolve!-Riccati in odenon1.
   symbolic if ode_order = 1 then
      ODESolve!-linear1(odecoeffs, driver, x)
   else or(
      ODESolve!-Run!-Hook(
         'ODESolve_Before_Lin_Hook,
         {odecoeffs,driver,y,x,ode_order,min_order}),
      ODESolve!-linearn
         (odecoeffs, driver, y, x, ode_order, min_order, nil),
      ODESolve!-Run!-Hook(
         'ODESolve_After_Lin_Hook,
         {odecoeffs,driver,y,x,ode_order,min_order}))$

algebraic procedure ODESolve!-linear!-basis!-recursive
   (odecoeffs, driver, y, x, ode_order, min_order);
   %% Always returns the solution in basis format.
   %% Internal linear solver called recursively.
   symbolic if ode_order = 1 then
      ODESolve!-linear1(odecoeffs, driver, x)
   else
      ODESolve!-linearn
         (odecoeffs, driver, y, x, ode_order, min_order, t)$


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Solve a linear first-order ODE by using an integrating factor.
% Based on procedure linear1 from module ode1ord by Malcolm MacCallum.

algebraic procedure ODESolve!-linear1(odecoeffs, driver, x);
   %% Solve the linear ODE reduced_ode = driver, where
   %% reduced_ode = A(x)*(dy/dx + P(x)*y), driver = A(x)*Q(x).
   %% Uses Odesolve!-Int to optionally turn off final integration.
   begin scalar A, P, Q;
      A := second odecoeffs;
      P := first odecoeffs/A;
      Q := driver/A;
      return if P then                  % dy/dx + P(x)*y = Q(x)
         begin scalar intfactor, !*combinelogs;
            traceode "It is solved by the integrating factor method.";
            %% intfactor simplifies better if logs are combined:
            symbolic(!*combinelogs := t);
            P := (P where tan(~x) => sin(x)/cos(x));
            intfactor := exp(int(P, x));
            return if Q then
               { {1/intfactor}, Odesolve!-Int(intfactor*Q,x)/intfactor }
            else { {1/intfactor} }
         end
      else <<                           % dy/dx = Q(x)
         traceode "It is solved by quadrature.";
         if Q then {{1}, Odesolve!-Int(Q, x)} else {{1}}
      >>
   end$


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Try to solve a linear ODE of order > 1.
% Based on procedure linearn from module linearn by Malcolm MacCallum.

% If the first integral of an exact ODE has constant coefficients, is
% of Euler type or has trivally reducible order then so does the
% original ODE.  Also, trivial order reduction preserves constant
% coefficients or Euler type, and the reduced ODE is not further
% reducible.  Hence, when the linear ODE solver is called recursively
% to solve a first integral of an exact ODE or a trivially order-
% reduced ODE, the argument `recursive' is used to avoid checking
% again whether it has constant coefficients, is of Euler type or has
% trivially reducible order.

algebraic procedure ODESolve!-linearn
   (odecoeffs, driver, y, x, ode_order, min_order, recursive);
   %% Solve the linear ODE: reduced_ode = driver.
   begin scalar lcoeff, odecoeffs1, driver1, solution;
      %% Make the ODE "monic" as assumed by some solvers:
      %% (Note that this makes algebraic factorization largely
      %% irrelevant!)
      if (lcoeff := part(odecoeffs, ode_order+1)) = 1 then <<
         odecoeffs1 := odecoeffs;
         driver1 := driver
      >> else <<
         odecoeffs1 := for each c in odecoeffs collect c/lcoeff; % low .. high
         %% Could discard last element of odecoeffs1 because it must be 1!
         driver1 := driver/lcoeff
      >>;
      if recursive then goto a;

      %% Test for constant coefficients:
      if odecoeffs1 freeof x then
         return ODESolve!-LCC(odecoeffs1, driver1, x, ode_order);

      traceode "It has non-constant coefficients.";

      %% Test for Euler form:
      if (solution :=
         ODESolve!-Euler(odecoeffs1, driver1, x, ode_order))
      then return solution;

      %% Test for trivial order reduction.  The result cannot have
      %% constant coeffs or Euler form, but it could be first order or
      %% exact or ...
      if min_order neq 0 and ode_order neq min_order and
         %% else would reduce to purely algebraic equation
         (solution := ODELin!-Reduce!-Order
            (odecoeffs, driver, y, x, ode_order, min_order))
      then return solution;

   a: %% Non-trivial solution techniques for recursive calls ...

      %% Test for exact form - try monic then original form:
      if (solution :=
         ODELin!-Exact(odecoeffs1, driver1, y, x, ode_order))
      then return solution;
      if lcoeff neq 1 and (solution :=
         ODELin!-Exact(odecoeffs, driver, y, x, ode_order))
      then return solution;

      %% Add other methods here ...
      %% Null return implies failure.

      %% FINALLY, test for a second-order special-function equation:
      if ode_order = 2 and
         (solution := ODESolve!-Specfn(odecoeffs1, driver1, x))
      then return solution
   end$


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Convert between basis and linear combination formats

% Solution basis output format:  { {B1, B2, ...}, PI }
% where {Bi} is a basis for the reduced ODE and PI is a particular
% intergral for the full ODE, which may be absent.
% This corresponds to the linear-combination output format
% y = ( for each B in {B1, B2, ...} sum newarbconst()*B ) + PI.

% This agrees with Maple, e.g. Maple V Release 5 gives
% > dsolve(diff(y(x),x) + y(x) = x, output=basis);
%                          [[exp(-x)], -1 + x]
% > dsolve(diff(y(x),x) + y(x) = x);
%                       y(x) = x - 1 + exp(-x) _C1

algebraic procedure ODESolve!-Basis2LinComb solution;
   %% Convert basis { {B1, B2, ...}, PI } to linear combination:
   begin scalar lincomb;
      lincomb := for each B in first solution sum <<newarbconst()>>*B;
      %% << >> above is NECESSARY to force immediate evaluation!
      if length solution > 1 then
         lincomb := lincomb + second solution;
      return lincomb
   end$

algebraic procedure ODESolve!-LinComb2Basis
      (lincomb, first_arb, last_arb);
   %% Convert linear combination to basis { {B1, B2, ...}, PI }:
   ODESolve!-LinComb2Basis1({}, lincomb, first_arb, last_arb)$

algebraic procedure ODESolve!-LinComb2Basis1
      (basis, lincomb, first_arb, last_arb);
   %% `basis' is a LIST of independent reduced_ode solutions.
   %% Algorithm is to recursively move components from lincomb to
   %% basis.
   begin scalar coeffs, C;  C := arbconst last_arb;
      coeffs := coeff(lincomb, C);
      if high_pow > 1 or smember(C, coeffs) then
         return                         % cannot convert
      else if high_pow = 1 then <<
         basis := second coeffs . basis;
         lincomb := first coeffs
      >>;
      %% else independent of this arbconst
      return if first_arb >= last_arb then
         { basis, lincomb }
      else
         ODESolve!-LinComb2Basis1
            (basis, lincomb, first_arb, last_arb-1)
   end$


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Solve a linear, constant-coefficient ODE.
% Based on code by Malcolm MacCallum.

% There are (at least) 4 ways to get the particular integral (P.I.)
% for a given driving term on the right of the equation:

% 1. The method of undetermined coefficients: this is similar to the
% integrator in that for a given driving term one has to find the
% functional form of the P.I. and then solve for the numerical
% coefficients in it.  Making it really general is as big a task as
% rewriting the integrator.

% 2. The method of variation of parameters: this expands the P.I. as a
% sum of functions of `x' times the linearly independent solutions in
% the complementary function (C.F.).

% 3. Factorise the linear operator (done anyway for the C.F.) and then
% apply for each root `m' the operation
%   ans := exp(m*x) * int(ans * exp(-m*x))
% This is a form of the "D-operator method".  N.B. Some `m' are
% complex.

% 4. Use Laplace transforms (and some kind of table lookup for the
% inverse transforms).

% The current implementation first tries to use the "D-operator
% method", but as soon as any integral fails to evaluate it switches
% to "variation of parameters".

algebraic procedure ODESolve!-LCC(odecoeffs, driver, x, ode_order);
   % Returns a solution basis or nil (if it fails).
   begin scalar auxvar, auxeqn, i, auxroots, solutions, PI;
      traceode "It has constant coefficients.";
      %% TEMPORARY HACK -- REBUILD AUXEQN:
      auxvar := symbolic gensym();  i := -1;
      auxeqn := for each c in odecoeffs sum c*auxvar^(i:=i+1);
      % First we solve for the auxiliary roots:
      auxroots := solve(auxeqn, auxvar);
      % and check the solution carefully:
      if ode_order neq
         (for each multi in multiplicities!* sum multi) then return
            traceode "But insufficient roots of auxiliary equation!";
      solutions := auxroots;
   a: if lhs first solutions neq auxvar then return
         traceode "But auxiliary equation could not be solved!";
      if (solutions := rest solutions) neq {} then goto a;
      % Now we find the complementary solution:
      solutions := ODESolve!-LCC!-CompSoln(auxroots, x);
      % Next the particular integral:
      if driver = 0 then return { solutions };
      if not (PI := ODESolve!-LCC!-PI(auxroots, driver, x)) then
         %% (Cannot use `or' as an algebraic operator!)
         PI := ODESolve!-PI(solutions, driver, x);
      return { solutions, PI }
   end$

algebraic procedure ODESolve!-LCC!-CompSoln(auxroots, x);
   %% Construct the complimentary solution (functions) from the roots
   %% of the auxiliary equation for a linear ODE with constant
   %% coefficients.  Pairs of complex conjugate roots are converted to
   %% real trigonometric form up to the minimum of their
   %% multiplicities (regardless of complex switch and parameters).
   %% `auxroots' is a list of equations with a temporary variable on
   %% the left and an auxilliary root on the right.  The root
   %% multiplicities are stored as a list in the global variable
   %% `multiplicities!*'.  `x' is the independent variable.
   begin scalar multilist, crootlist, ans, multi, imroot, exppart;
      %% crootlist will be a list of lists of the form
      %% {unpaired_complex_root, multiplicity}.
      multilist := multiplicities!*;
      crootlist := {};  ans := {};
      for each root in auxroots do <<
         root := rhs root;
         multi := first multilist;  multilist := rest multilist;
         % Test for complex roots:
         imroot := impart!* root;
         if imroot = 0 then <<
            exppart := exp(root*x);
            for j := 1 : multi do ans := (x**(j-1)*exppart) . ans
         >> else
         begin scalar conjroot, conjmulti;
            %% Cannot assume anything about the order of the roots in
            %% auxroots, so build a list of the complex roots found to
            %% avoid using complex conjugate pairs twice.
            conjroot := conj!* root;  conjmulti := 0;
            %% Essentially do assoc followed by delete if found:
            crootlist := for each root in crootlist join
               if first root = conjroot then <<
                  conjmulti := second root; {}
               >> else {root};
            if conjmulti then           % conjugate pair found:
            begin scalar minmulti;
               exppart := exp(repart!* root*x);
               minmulti := min(multi, conjmulti);
               imroot := abs imroot;    % to avoid spurious minus sign
               imroot := (imroot where abs ~x => x); % to avoid spurious abs!
               for j := 1 : minmulti do
                  ans := (x**(j-1)*cos(imroot*x)*exppart) .
                     (x**(j-1)*sin(imroot*x)*exppart) . ans;
               if multi neq conjmulti then <<
                  %% Skip this unlikely case if possible
                  minmulti := minmulti + 1;
                  exppart := exp(root*x);
                  for j := minmulti : multi do
                     ans := (x**(j-1)*exppart) . ans;
                  exppart := exp(conjroot*x);
                  for j := minmulti : conjmulti do
                     ans := (x**(j-1)*exppart) . ans
               >>
            end
            else crootlist := {root, multi} . crootlist
         end
      >>;
      %% Finally include unpaired complex roots:
      for each root in crootlist do <<
         exppart := exp(first root*x);
         multi := second root;
         for j := 1 : multi do ans := (x**(j-1)*exppart) . ans
      >>;
      return ans
   end$

% The following procedures process complex-valued expressions with
% regard only to their EXPLICIT complexity, i.e. assuming that all
% symbolic quantities are pure real.  They need to work with the
% complex switch both on and off, which is slightly tricky!

algebraic(vars!-are!-real := {repart ~x => x, impart ~x => 0})$

algebraic procedure repart!* u;
   << u := repart u;  u where vars!-are!-real >>$

algebraic procedure impart!* u;
   << u := impart u;  u where vars!-are!-real >>$

algebraic procedure conj!* u;
   << u := conj u;  u where vars!-are!-real >>$

algebraic procedure ODESolve!-LCC!-PI(auxroots, driver, x);
   % Try to construct a particular integral using the `D-operator
   % method'.  Factorise the linear operator (done anyway for the
   % C.F.) and then apply for each root m the operation
   % ans := exp(m*x) * int(ans * exp(-m*x));
   % N.B. Some m may be complex.
   % See e.g. Stephenson, section 21.8, p 410.
   % Returns nil if any integral cannot be evaluated.
   begin scalar exp_mx, multiplicities, multi;
      traceode
         "Constructing particular integral using `D-operator method'.";
      multiplicities := multiplicities!*;
      while driver and auxroots neq {} do <<
         exp_mx := exp((rhs first auxroots)*x);
         driver := driver/exp_mx;
         multi := first multiplicities;
         while driver and multi >= 1 do <<
            driver := int(driver, x);
            if driver freeof int then multi := multi - 1 else driver := 0
         >>;
         driver := exp_mx*driver;
         auxroots := rest auxroots;
         multiplicities := rest multiplicities
      >>;
      if driver = 0 then
         traceode "But cannot evaluate the integrals, so ..."
      else return driver
   end$

algebraic procedure ODESolve!-PI(solutions, R, x);
   % Given a "monic" forced linear nth-order ODE

   %   y^(n) + a_(n-1)(x)y^(n-1) + ... + a_1(x)y = R(x)

   % and a set of n linearly independent solutions {yi(x)} to the
   % unforced ODE, construct a particular solution of the forced ODE
   % in the form of a (single) integral representation by the method
   % of variation of parameters.
   begin scalar n;
      traceode
         "Constructing particular integral using `variation of parameters'.";
      return
         if (n := length solutions) = 2 then
         begin scalar y1, y2, W;
            y1 := first solutions;  y2 := second solutions;
            %% The Wronskian, kept separate to facilitate tracing:
            W := trigsimp(y1*df(y2, x) - y2*df(y1, x));
            traceode "The Wronskian is ", W;
            R := R/W;
            return -ode!-int(y2*R, x)*y1 + ode!-int(y1*R, x)*y2
         end
         else
         begin scalar Wmat, ys, W, i;
            %% Construct the (square) Wronskian matrix of the solutions:
            Wmat := {ys := solutions};
            for i := 2 : n do
               Wmat := (ys := for each y in ys collect df(y,x)) . Wmat;
            load_package matrix;              % to define mat
            Wmat := list2mat reverse Wmat;
            %% The Wronskian (determinant), kept separate for tracing:
            W := trigsimp det Wmat;
            traceode "The Wronskian is ", W;
            R := R/W;  i := 0;
            return
               for each y in solutions sum
                  ode!-int(cofactor(Wmat, n, i:=i+1)*R, x) * y
         end
   end$

% This facility should be in the standard matrix package!
symbolic operator list2mat$
symbolic procedure list2mat M;
   % Input:  (list (list A B ...) (list C D ...) ...)
   % Output: (mat  (A B ...) (C D ...) ...)
   'mat . for each row in cdr M collect cdr row$


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Special cases of non-constant coefficients:

algebraic procedure ODESolve!-Euler(odecoeffs, driver, x, ode_order);
   %% Solve a (MONIC) ODE having (essentially) the form
   %% reduced_ode = x^n df(y,x,n) + ... + a_{n-1} x df(y,x) + a_n y = driver
   %% odecoeffs = {a_n, a_{n-1} x, ..., a_0 x^n} / (a_0 x^n)
   %%           = {a_n/(a_0 x^n), a_{n-1}/(a_0 x^{n-1}), ..., a_1/(a_0 x), 1}
   begin scalar tmp, shift, i, c, solution;
      odecoeffs := reverse odecoeffs;   % high .. low

      %% Check for possible Euler or "shifted Euler" form.
      %% Find second non-zero ode coeff:
      tmp := rest odecoeffs;  i := 1;
      while first tmp = 0 do << tmp := rest tmp; i := i+1 >>;
      tmp := first tmp;                 % second non-zero ode coeff
      tmp := den tmp;                   % ax^i or a(x+b)^i
      tmp := reverse coeff(tmp, x);     % high .. low
      if high_pow neq i then return;    % not Euler
      if second tmp then <<             % "shifted Euler"
         shift := second tmp/(i*first tmp); % b
         driver := sub(x=x-shift, driver) % x -> x-b
      >>;

      tmp := {first odecoeffs};  i := 0;
      odecoeffs := rest odecoeffs;
   a: if odecoeffs neq {} then <<
         c := first odecoeffs * (x+shift)^(i:=i+1);
         if not(c freeof x) then return; % not Euler
         tmp := c . tmp;
         odecoeffs := rest odecoeffs;
         go to a
      >>;
      odecoeffs := tmp;

      traceode "It is of the homogeneous (Euler) type ",
         if shift then "(with shifted coefficients) " else "",
         "and is reducible to a simpler ODE ...";
      i := -2;
      tmp := for each c in odecoeffs sum <<
         i := i + 1;
         c * for j := 0 : i product (x-j)
      >>;
      odecoeffs := coeff(tmp, x);       % TEMPORARY HACK!
      driver := sub(x=e^x, driver*x^ode_order);
      solution := ODESolve!-LCC(odecoeffs, driver, x, ode_order);
      solution := sub(x=log x, solution);
      if shift then solution := sub(x=x+shift, solution);
      return solution
   end$

algebraic procedure ODELin!-Exact(P_list, driver, y, x, n);
   %% Computes a (linear) first integral if ODE is an exact linear
   %% n'th order ODE P_n(x) df(y,x,n) + ... + P_0(x) y = R(x).
   begin scalar P_0, C, Q_list, Q, const, soln, PI;
      P_0 := first P_list;
      P_list := reverse rest P_list;         % P_n, ..., P_1
      %% ODE is exact if C = df(P_n,x,n) - df(P_{n-1},x,{n-1}) + ...
      %% + (-1)^{n-1} df(P_1,x) + (-1)^n P_0 = 0.
      for each P in P_list do C := P - df(C,x);
      C := P_0 - df(C,x);               % C = 0 if exact
      if C then return;
      Q_list := {};
      for each P in P_list do
         Q_list := (Q := P - df(Q,x)) . Q_list; % Q_0, ..., Q_{n-1}
      driver := int(driver, x) + (const := symbolic gensym());
      %% The first integral is the LINEAR (n-1)'th order ODE
      %% Q_{n-1}(x) df(y,x,n) + ... + Q_0(x) y = int(R(x),x).
      traceode "It is exact, and the following linear ODE of order ",
         n-1, " is a first integral:";
      if symbolic !*trode then <<
         C := y;
         soln := first Q_list*y +
            ( for each Q in rest Q_list sum Q*(C := df(C,x)) );
         write soln = driver
      >>;
      %% Recurse on the order:
      C := Q_list;
      %% First-integral ODE must have min order 0, since input ODE was
      %% already order-reduced.
      soln := ODESolve!-linear!-basis!-recursive
         (Q_list, driver, y, x, n-1, 0);
      PI := second soln;                % MUST exist since driver neq 0
      PI := coeff(PI, const);           % { real PI, extra basis fn }
      return if high_pow = 1 then
         if first PI then { second PI . first soln, first PI }
         else { second PI . first soln }
      else <<
         %% This error should now be redundant!
         write "*** Internal error in ODELin!-Exact:",
            " cannot separate basis functions! ";
         write "(Probably caused by `noint' option.)";
         soln
      >>
   end$

algebraic procedure ODELin!-Reduce!-Order
   (odecoeffs, driver, y, x, ode_order, min_order);
   %% If ODE does not explicitly involve y (and perhaps low order
   %% derivatives) then simplify by reducing the effective order
   %% (unless there is only one) and try to solve the reduced ODE to
   %% give a first integral.  Applies only to ODEs of order > 1.
   begin scalar solution, PI;
      ode_order := ode_order - min_order;
      for ord := 1 : min_order do odecoeffs := rest odecoeffs;
      traceode "Performing trivial order reduction to give the order ",
         ode_order, " linear ODE with coefficients (low -- high): ",
         odecoeffs;
      solution := ODESolve!-linear!-basis!-recursive
         (odecoeffs, driver, y, x, ode_order, 0);
      if not solution then <<
         traceode "But ODESolve cannot solve the reduced ODE! ";
         return
      >>;
      traceode "Solution of order-reduced ODE is ", solution;
      traceode "Restoring order, ", y => df(y,x,min_order),
         ", to give: ", df(y,x,min_order) = solution,
         " and re-solving ...";
      %% = lin comb of fns of x, so just integrate min_order times:
      if length solution > 1 then       % PI = particular integral
         PI := second solution;
      solution := append(
         for each c in first solution collect
            ODESolve!-multi!-int(c, x, min_order),
         %% and add min_order extra basis functions:
         for i := min_order-1 step -1 until 0 collect x^i );
      return if PI then
         { solution, ODESolve!-multi!-int(PI, x, min_order) }
      else
         { solution }
   end$

endmodule$

end$