Perturbation

Stephan Kaufmann

Perturbation expansions for the solutions of ordinary differential equations.

Reference

Author

Stephan Kaufmann
Institute of Mechanics
Swiss Federal Institute of Technology
ETH Zentrum, HG F38.4 
CH-8092 Zuerich
Switzerland
kaufmann@ifm.ethz.ch

Copyright Notice

Copyright © 1992 by Stephan Kaufmann

Permission is hereby granted to make copies of	this file for any purpose other than direct profit,	or as part of a commercial product, provided this copyright notice is left intact.	Sale is prohibited.
	
Permission is hereby granted to reproduce part or all of this file provided that the source is acknowledged.

Version

Notebook Version 2.2 (E)

History

Changes in Version 2.2 (August 93):
Automatic simplification added (for Mma 2.2)

Requirements

This package is based on SKPackages`NonNegativeQ`. This package must be present.

Mathematica Version: 2.0 or higher.

If[$VersionNumber < 2.0, Print[
  "This package requires Mathematica version 2.0 or higher."]]

Warning

The function Limit is extended by a further rule.

Assumptions

In this package, all symbols are assumed to be real.

In the function MatchedSolutionList, the small parameter eps and the independent variables t and tau are assumed to be positive (and real).

Discussion

Literature

Main Literature

Cole, J. D.: Perturbation Methods. Blaidell, Waltham, Mass.
Nayfeh, A. H.: Perturbation Methods. Wiley, New York, 1973.
Nayfeh, A. H.: Introduction to Perturbation Techniques. Wiley, New York, 1981.

Further Literature

Kirchgraber, U. and E. Stiefel: Methoden der analytischen Störungsrechnung und ihre Anwendungen. Teubner, Stuttgart, 1978.
Rand, H. R. and D. Armbruster: Perturbation Methods, Bifurcation Theory and Computer Algebra. Springer 1987.

Overview

The aim of this package is to show a possible implementation of perturbation methods with Mathematica. It can be used to generate educational examples of perturbation expansions.

The current version of the package covers asymptotic expansions of the solutions of ordinary differential equations (ODE's) with respect to polynomial gauge functions (according to the main sources of literature cited above). The methods of straight forward expansions, strained coordinates and matched and composite solutions are implemented.

The method of strained coordinates only uses the simple scheme of Linstedt-Poincaré.

The matching scheme used for matched and composite expansions equals the n-th order asymptotic expansion of the inner solution of order n (transformed to outer coordinates) to the n-th order asymptotic expansion of the outer solution of order n (transformed to inner coordinates).

Examples

Auxiliary Functions

The function PolyOrderList[equations, x[t], eps, maxOrder] substitutes Sum[eps^i x[i][t], {i, 0, maxOrder}] for x[t] into the differential equation and the initial conditions, and calculates the resulting equations of orders 0 to maxOrder.

PolyOrderList[{x''[t] + eps (x'[t])^3 + x[t] == 0,
  x[0] == a, x'[0] == b}, x[t], eps, 3]

It is sometimes necessary to split a combined list of differential equations and initial conditions into those two parts:

Separate[{x''[t] + eps (x'[t])^3 + x[t] == 0,
  x[0] == a, x'[0] == b}, x[t], t]

The following function returns the order of the highest derivative in a list of equations:

MaxDerivative[{x''[t] + eps (x'[t])^3 + x[t] == 0,
  x[0] == a, x'[0] == b}]

Asymptotic Expansions

The coefficient c[i] of the asymptotic expansion of a function f[eps] (for eps -> 0) with respect to a list g of gauge functions (in eps) is calculated (for n = 1, ... , m) according to:

c[n] = Limit[(f[eps] - Sum[c[i] g[[i]], {i, n-1}])/g[[n]],
  eps -> 0, Direction -> -1].

The function AsCoefficients[f[eps, ...], gaugeFcts, eps] returns the list of the coefficients {c[1], ... , c[m]}:

AsCoefficients[Sin[eps + t], Table[eps^j, {j, 0, 6}], eps]

AsExpansion[f[eps, ...], gaugeFcts, eps] calculates the asymptotic expansion of the given function f, i.e.: Sum[c[i] g[[i]], {i, m}].

AsExpansion[Sin[eps + x], Table[eps^j, {j, 0, 6}], eps]

Regular Perturbation, Straight Forward Expansions

The asymptotic expansion of the solution of ordinary differential equations is generated by substituting Sum[eps^i x[i][t], {i, 0, maxOrder}] for x[t] into the equation and solving the resulting equations (in increasing order).

PolyOrderSolutionList[equations, x[t], t, eps, maxOrder] returns the list of the asymptotic expansions of the solution (of the equations equations, in the function x[t], with the small parameter eps) up to the desired maximum order. Initial conditions are passed as equations, analogously to the parameters of DSolve.

PolyOrderSolutionList[
  { x''[t] + eps (x'[t])^2 + x[t] == 0,
    x[0] == 1, x'[0] == 0 },
  x[t], t, eps, 2]

Strained Coordinates

The method of strained coordinates is useful when the exact solution has a different frequency than the solution of the unperturbed system (eps = 0). Straightforward expansions of such equations produce secular terms (growing polynomially in the independent variable) and therefore are not uniform asymptotic expansions for infinite domains.

The method of Linstedt-Poincaré works with an almost identical transformation of the independent variable t to an new variable s, according to:

t = s (1 + Sum[eps^i o[i], {i, 0, maxOrder}]).

The o[i] are chosen such that the secular terms vanish.

The following straight forward expansion has such a secular term (6 t eps Sin[t]):

PolyOrderSolutionList[
  {x''[t] + x[t] + eps x[t]^3 == 0, x[0] == 1, x'[0] == 0}, 
    x[t], t, eps, 1]

The method of strained coordinates results in a better approximation. StrainedSolutionList[equations, x[t], t, eps, maxOrder] returns this expansion (the parameters are analogous to PolyOrderSolutionList).

StrainedSolutionList[
  {x''[t] + x[t] + eps x[t]^3 == 0, x[0] == 1, x'[0] == 0}, 
    x[t], t, eps, 1]

Matched and Composite Expansions

If the small parameter (eps) multiplies the highest derivative, a straightforward expansion does not work, because there are not enough free constants to satisfy the initial conditions. Physically spoken, there is a boundary layer, where the changes of the functions have a different scale. The method of matched and composite expansions enlarges this boundary layer with a suitable coordinate transformation (e.g.: tau = t/eps). It solves the equations in the bounary layer and outside it by asymptotic expansions. The solutions are matched such that they agree in an overlapping region, and then combined.

Such expansions can display the physical properties of solutions more clearly than the formulas for the exact solution (if available).

Of course, the asymptotic expansion of the exact solution should be asymptotically equivalent to the composite solution.

The following example shows a simple equation which can not be expanded in a straightforward fashion:

PolyOrderSolutionList[{eps x''[t] + x'[t] + x[t] == 0, 
  x[0] == a, x[1] == b}, x[t], t, eps, 1]

The auxiliary function BoundaryLayerEquation[equation, bndRel, x[t], t, tau, eps] calculates the equation for the boundary layer by transforming the independent variable t to tau, where t and tau are connected according to a given equation, e.g.: tau == t/eps. For a small eps, this enlarges the interval from t=0 to a small positive value of t.

BoundaryLayerEquation[eps x''[t] + x'[t] + x[t]  == 0, 
  tau == t/eps, x[t], t, tau, eps]

The function MatchedSolutionList[equations, bndRel, x[t], t, tau, eps, maxOrder] returns the list of the solutions up to the desired order.

MatchedSolutionList[{eps x''[t] + x'[t] + x[t] == 0, 
  x[0] == a, x[1] == b},
  tau == t/eps, x[t], t, tau, eps, 2]

We can choose different initial conditions:

MatchedSolutionList[{eps x''[t] + x'[t] + x[t] == 0, 
  x[0] == 1, x'[0] == 0},
  tau == t/eps, x[t], t, tau, eps, 2]

In this following example, DSolve gives us the exact solution:

DSolve[{eps x''[t] + x'[t] + x[t] == 0, 
  x[0] == 1, x'[0] == 0}, x[t], t]

Notice that the dependence on the small parameter eps can be seen much more clearly in the expansion!

Implementation

Preparation

BeginPackage["SKPackages`Perturbation`",
  "SKPackages`NonNegativeQ`"];

Separate::usage =
"Separate[equations, x[t], t] separates the given list of differential equations (in the function x[t] with the independent variable t) and initial conditions into a list of differential equations and a list of initial conditions.";

MaxDerivative::usage =
"MaxDerivative[equ] returns the highest derivative appearing in the given equation or list of equations.";

AsCoefficients::usage =
"AsCoefficients[f[eps, ...], gaugeFcts, eps] calculates the asymptotic expansion of the function f[eps, ...] with respect to the list of gauge functions gaugeFcts and returns the list of coefficients corresponding to the list of gauge functions. The gauge functions are functions of the small positive parameter eps. The expansion is asymptotic in the limit eps -> 0.\n
Example:\n
AsCoefficients[Sin[eps + t], Table[eps^j, {j, 0, 6}], eps]";

AsExpansion::usage =
"AsExpansion[f[eps, ...], gaugeFcts, eps] calculates the asymptotic expansion of the function f[eps, ...] with respect to the list of gauge functions gaugeFcts. The gauge functions are functions of the small positive parameter eps. The expansion is asymptotic in the limit eps -> 0.\n
Example:\n
AsExpansion[Sin[eps + x], Table[eps^j, {j, 0, 6}], eps]";

PolyOrderList::usage =
"PolyOrderList[equationList, x[t], eps, maxOrder] substitutes Sum[eps^i x[i][t0], {i, 0, maxOrder}] in to the list of ordinary differential equations and initial conditions equationList. It extracts the list equations of rising orders of eps up to maxOrder.\n
Example:\n
PolyOrderList[{x''[t] + x[t] + eps x[t]^3 == 0, x[0] == 1, x'[0] == eps}, x[t], eps, 2]";

PolyOrderSolutionList::usage =
"PolyOrderSolutionList[equations, x[t], t, eps, maxOrder] calculates an asymptotic expansion of the solution of the given ordinary differential equation(s) containing the functions x[t] of the independent variable t and a small parameter eps. The solution is generated by substituting Sum[eps^i x[i][t], {i, 0, maxOrder}] into the equation(s) and solving for increasing orders in eps. The list of solutions up to the chosen order maxOrder is calculated.\n
Example:\n
PolyOrderSolutionList[{x''[t] + eps (x'[t])^3 + x[t] == 0, x[0] == 1, x'[0] == 0}, x[t], t, eps, 2]";

StrainedSolutionList::usage =
"StrainedSolutionList[equations, x[t], t, eps, maxOrder] calculates an asymptotic expansion of the solution of the given ordinary differential equation(s) containing the functions x[t] of the independent variable t and a small parameter eps. It uses the method of strained coordinates of Linstedt-Poincare. The solution is generated by substituting s*(Sum eps^i o[i], {i, maxOrder}) for the time t into the equations and Sum[eps^i x[i][s], {i, 0, maxOrder}] into the resulting equation(s). The o[i] are chosen such that the secular terms vanish. The list of solutions up to the chosen order maxOrder is calculated.\n
Example:\n
StrainedSolutionList[
  {x''[t] + x[t]+ eps x[t]^3 == 0, x[0] == 1, x'[0] == 0}, x[t], t, eps, 1]";

BoundaryLayerEquation::usage =
"BoundaryLayerEquation[equation, bndRel, x[t], t, tau, eps] calculates the boundary layer equation corresponding to the given equation, which contains the function x[t] in the variable t. The boundary variable is called tau. The relation between t and tau is given as an equation. The small parameter is called eps.\n
Example:\n
BoundaryLayerEquation[eps x''[t] + x'[t] + x[t]  == 0, tau == t/eps, x[t], t, tau, eps]";

MatchedSolutionList::usage =
"MatchedSolutionList[equations, bndRel, x[t], t, tau, eps, maxOrder] calculates the asymptotic expansion of the solution of the given ordinary differential equation (with initial conditions), containing the functions x[t] of the independent variable t and the small parameter eps. The method of matched composite solutions is applied. The boundary layer equation is generated by substituting the boundary variable tau for t into the equation. The relations between t and tau ist given by the equation bndRel. The equation in the boundary layer (inner equation) and the given equation (outer equation) are solved and matched by the setting the innner (asymptotic) expansion of the outer solution equal to the outer (asymptotic) expansion of the inner solution, for all orders. The result ist returned as a list of solutions up to the given order maxOrder.\n
Example:\n
MatchedSolutionList[{eps x''[t] + x'[t] + x[t] == Sin[t], x[0] == 0, x'[0] == 1}, tau == t/eps, x[t], t, tau, eps, 1]";

Perturbation::inising =
"The differential equation of order `1` does not have a solution which satisfies the given initial conditions. Consider attempting a singular perturbation!";

Perturbation::noinner =
"The inner solution of order `1` could not be found.";

Perturbation::noouter =
"The outer solution of order `1` could not be found.";

Perturbation::nomatch =
"Matching could not be achieved with the given boundary variable.";

Perturbation::rec =
"The implemented method of strained coordinates leads to an infinite recursion.";

Asymptotic::notfound =
"A limit could not be determined. You can try to extend the Limit function with suitable definitions.";

Begin["`Private`"];

Definitions

Extension of Limit

prot = Unprotect[Limit];

Limit[a_. Exp[b_. x_^n_] + c_., x_ -> 0,
    Direction -> d_] := 
  Limit[c, x -> 0, Direction -> d] /; 
    NonNegativeQ[-n] && FreeQ[b, x] && 
    (!OddQ[n] && NonNegativeQ[-b] ||
     OddQ[n] && NonNegativeQ[d b]) &&
    ( FreeQ[a, x] ||
      PolynomialQ[Numerator[a], x] &&
        PolynomialQ[Denominator[a], x] ||
      Limit[a, x -> 0] < Infinity)

Protect[Evaluate[prot]];

Auxiliary Functions

Separate[equations_List, x_, t_] :=
  {Select[equations, !FreeQ[#, t]&], Select[equations, FreeQ[#, t]&]}

MaxDerivative[equ_] := 
  Max[Cases[equ, Derivative[_][_][_], Infinity] /. 
    Derivative[n_][_][_] -> n] /. -Infinity -> 0;

PolyOrderList[equations_List, x_, eps_, maxOrder_Integer] :=
  Module[{j, h = Head[x], substEqu},
    substEqu = equations /. 
      {
        h[t_] -> Sum[eps^i h[i][t], {i, 0, maxOrder}],
        Derivative[n_][h][t_] ->
          Sum[eps^i Derivative[n][h[i]][t], {i, 0, maxOrder}]
      } // ExpandAll;
    Table[Map[Coefficient[#, eps, j]&, substEqu, {2}], {j, 0, maxOrder}]
  ]

Asymptotic Expansions

AsCoefficients[f_, gaugeFcts_List, eps_] :=
  Module[
    { c, i = 1, 
      tempSum = 0, lastIndex = Length[gaugeFcts]},
    While[i < lastIndex,  
      c[i] = Limit[Expand[(f - tempSum)/gaugeFcts[[i]]], 
        eps -> 0, Direction -> -1];
      If[!FreeQ[c[i], Limit],
        Message[Asymptotic::notfound]; Return["error"]];
      tempSum += c[i] gaugeFcts[[i]];
      i++
    ];
    c[lastIndex] = Limit[Expand[
      (f - tempSum)/gaugeFcts[[lastIndex]]],
      eps -> 0, Direction -> -1];
    If[!FreeQ[c[i], Limit],
      Message[Asymptotic::notfound]; Return["error"],
      Return[Array[c, lastIndex]]
    ]
  ]

AsExpansion[f_, gaugeFcts_List, eps_] :=
  Module[{coeffs = AsCoefficients[f, gaugeFcts, eps]},
    If[!FreeQ[coeffs, "error"],
      Return["error"],
      Return[coeffs . gaugeFcts]
    ]
  ]

Regular Perturbation, Straight Forward Expansions

PolyOrderSolutionRuleListPure[equations_List, x_, t_, eps_, 
    maxOrder_Integer] :=
  Module[{h = Head[x], SolveOrder},
    SolveOrder[rule_, equation_] :=
      Module[{i, res},
        If[!FreeQ[rule, "error"], Return["error"]];
        i = Length[rule];
        Print["Solving equation of order ", i, " ..."];
        res = DSolve[equation /. rule, h[i], t] //
          Simplify;
        Which[
          !FreeQ[res, DSolve], 
            Return["error"],
          res == {},
            Message[Perturbation::inising, i];
            Return["error"],
          True,
            Append[rule, res[[1, 1]]]
        ]
      ];
      Fold[SolveOrder, {},
        PolyOrderList[equations, x, eps, maxOrder]]
  ]

PolyOrderSolutionList[equations_List, x_, t_, eps_, 
    maxOrder_Integer] :=
  Module[
    {
      h = Head[x],
      pure = PolyOrderSolutionRuleListPure[equations, 
        x, t, eps, maxOrder]
    },
    If[!FreeQ[pure, "error"], Return["error"]];
    Drop[FoldList[Plus, 0, 
      Table[eps^i h[i][t] /. pure, {i, 0, maxOrder}]], 1]
  ]

Strained Coordinates

o::usage = "Private variable";

StrainedSolutionRuleListPure[equations_List, x_, t_, eps_, 
    maxOrder_Integer] :=
  Module[
    {
      SolveStrainedOrder, MaxDer, substEqu, s, om,
      tmpRes, invOm,
      h = Head[x],
      i = 0,      
    },
    SolveStrainedOrder[rules_, equation_] :=
      Module[{i, diff, sub},
        If[!FreeQ[rules, "error"], Return["error"]];
        i = Length[rules[[1]]];
        Print["Solving equation of order ", i, " ..."];
        diff = DSolve[equation /. Flatten[rules], h[i], s] //
          Simplify;
        Print["Removing secular terms of order ", i, " ..."];
        If[FreeQ[diff, DSolve],
          sub = Solve[Coefficient[
            h[i][s] /. diff[[1]] // ExpandAll, s] == 0,
            o[i]],
          Return["error"]
        ];
        Which[
          diff === {},
            Message[Perturbation::inising, i];
            Return["error"],
          !FreeQ[sub, s],
            Message[Perturbation::rec];
            Return["error"],
          sub === {} || sub === {{}}, 
            {Append[rules[[1]], diff[[1, 1]]], rules[[2]]},
          True,
            MapThread[Append,
              {rules, {diff[[1, 1]] /. sub[[1, 1]], sub[[1, 1]]}}]
        ]
      ];
            
    om = 1 + Sum[eps^j o[j], {j, maxOrder}];
    
    (* substitute t -> s (1 + Sum[eps^j o[j], {j, maxOrder}]) *)
    substEqu = Map[Together, equations /. 
      {
        Derivative[n_][h][t] -> om^-n Derivative[n][h][s],
        Derivative[n_][h][t0_] ->
          ((om^-n Derivative[n][h][s])/. s -> t0/om),
        h[t] -> h[s],
        h[t0_] -> h[t0/om],
        t -> s om
      },
      {2}];
      
    (* multiply with denominators *)
    substEqu = Thread[
      #(1 + Sum[eps^j o[j], {j, maxOrder}])^MaxDerivative[#],
        Equal]& /@ substEqu;
      
    (* calculate {rules for x[i], rules for o[i]}: *)
    tmpRes = Fold[SolveStrainedOrder, {{}, {}}, 
      PolyOrderList[substEqu, x, eps, maxOrder]];
      
    (* catch errors *)
    If[!FreeQ[tmpRes, "error"], Return["error"]];
      
    (* expand 1/om *)
    invOm = ExpandAll[Normal[Series[1/om,
      {eps, 0, maxOrder}]]];
    
    (* generate list of incr. order in eps of expansion of t/om *)
    invOm = t FoldList[(#1 + #2 Coefficient[invOm, #2])&, 1, 
      Table[eps^i, {i, maxOrder}]];
      
    (* return {rules for x[i], orders of t/om} *)
    Return[{tmpRes[[1]], invOm /. tmpRes[[2]]}]
  ]

StrainedSolutionList[equations_List, x_, t_, eps_, 
    maxOrder_Integer] :=
  Module[
    {
      h = Head[x],
      res,
      pure = StrainedSolutionRuleListPure[equations, 
        x, t, eps, maxOrder]
    },
    If[!FreeQ[pure, "error"], Return["error"]];
    res = Drop[FoldList[Plus, 0, 
      Table[eps^i h[i][pure[[2, i+1]]] /.
        pure[[1]], {i, 0, maxOrder}]], 1];
    If[!FreeQ[res, o],
      Message[Perturbation::rec]; Return["error"],
      Return[res]]
  ]

Matched and Composite Expansions

SetAttributes[BoundaryLayerEquation, Listable]

BoundaryLayerEquation[equation_Equal, bndRel_Equal,
  x_, t_, tau_, eps_] :=
  Module[
    { 
      res,
      h = Head[x],
      rel = tau /. Solve[bndRel, tau][[1]]
    },
    res = equation /.
      { 
        x -> h[tau],
        Derivative[n_][h][t] -> D[rel, t]^n Derivative[n][h][tau]
      } /. Solve[bndRel, t][[1]];
    Thread[res Denominator[Together[res[[1]] - res[[2]]]], Equal] //
      ExpandAll
  ]

MatchedSolutionList[equations_List, bndRel_Equal,
    x_, t_, tau_, eps_, maxOrder_] :=
  Module[
    {
      gaugeList, constList, CO, CI,
      h = Head[x],
      locT, locTau, locEps, setLocalRule, setGlobalRule,
      innerEqus, outerEqus, inits, subInits,
      SolveMatchedOrder, pureSolutionList, ieo, solutionList
    },
    
    Print["Preparing ..."];
    Off[Solve::svars];
    
    (* List of gauge Functions *)
    gaugeList = Table[locEps^n, {n, 0, maxOrder}];
    
    (* Integration constants for inner and outer equation *)
    constList = Flatten[Table[
      {CO[i, j], CI[i, j]}, {i, 0, maxOrder},
        {j, MaxDerivative[equations]}]];
        
    (* Set (locally substituted) independent variables and eps positive *)
    Positive[locT] ^= True;
    Positive[locTau] ^= True;
    Positive[locEps] ^= True;
    setLocalRule = {t -> locT, tau -> locTau, eps -> locEps};
    setGlobalRule = {locT -> t, locEps -> eps};
    
    (* Calculate outer and inner equation and initials *)
    outerEqus = PolyOrderList[Separate[equations, h[t], t][[1]],
      h[t], eps, maxOrder] /. setLocalRule;
    innerEqus = PolyOrderList[BoundaryLayerEquation[
      Separate[equations, h[t], t][[1]], bndRel, h[t], t, tau, eps],
      h[t], eps, maxOrder] /. setLocalRule;
    inits = Separate[equations, h[t], t][[2]];
    
    (* Calculate solution and matching condition of one order *)
    SolveMatchedOrder[rules_, ioequations_] :=
      Module[{i, outerRes, innerRes, match},
        
        If[!FreeQ[rules, "error"], Return["error"]];
        i = Length[rules[[1]]];
          
        Print["Solving outer equation of order ",
          i, " ..."];
        outerRes = (DSolve[ioequations[[1]] /. rules[[1]], h[i], locT] /.
          C[k_] -> CO[i, k]) // Simplify;
        If[!FreeQ[outerRes, DSolve],
          Message[Perturbation::noouter, i];
          Return["error"]];
          
        Print["Solving inner equation of order ", i, " ..."];
        innerRes = (DSolve[ioequations[[2]] /.
          rules[[2]], h[i], locTau] /.
          C[k_] -> CI[i, k]) // Simplify;
        If[!FreeQ[innerRes, DSolve],
          Message[Perturbation::noinner, i];
          Return["error"]];
          
        Print["Matching order ", i, " ..."];
        match = Solve[
          AsCoefficients[locEps^i h[i][locTau locEps] /.
              outerRes[[1]], 
            Take[gaugeList, i+1], locEps][[i+1]] ==
          AsCoefficients[locEps^i h[i][locT/locEps] /.
              innerRes[[1]], 
            Take[gaugeList, i+1], locEps][[i+1]],
          constList];          
        If[match == {{}} || !FreeQ[match, Limit] || 
            !FreeQ[match, Solve],
          Message[Perturbation::nomatch];
          Return["error"]
        ];

        MapThread[Append, {rules, 
          { outerRes[[1, 1]], innerRes[[1, 1]], 
            match[[1, 1]]}}
        ]
      ];
      
      (* List of pure inner and outer solutions and matching rules *)
      pureSolutionList = Fold[SolveMatchedOrder, 
        {{}, {}, {}}, Transpose[{outerEqus, innerEqus}]];
      
      (* catch errors *)
      If[!FreeQ[pureSolutionList, "error"],
        Return["error"]];

      Print["Calculating composite solution ..."];
      
      (* Inner expansion of outer solution *)
      ieo =
        AsCoefficients[Sum[locEps^i h[i][locT], {i, 0, maxOrder}] /. 
          pureSolutionList[[1]] /. locT -> 
            locTau locEps, gaugeList, locEps
        ] /. locTau -> locT/locEps;
        
      (* Composite solution *)
      solutionList = 
        Drop[FoldList[
          Plus, 0, gaugeList*
            ( (Table[h[i][locT], {i, 0, maxOrder}] /.
                  pureSolutionList[[1]]) +
              (Table[h[i][locTau /. Solve[bndRel /.
                  setLocalRule, locTau][[1]]] /. 
                  pureSolutionList[[2]], {i, 0, maxOrder}]) - 
              ieo) /.
            pureSolutionList[[3]] // ExpandAll
          ], 1];
          
      Print["Setting initials ..."];
      
      (* Substitute solution into initials *)
      subInits = ((Thread[# locEps^MaxDerivative[#], 
          Equal])& /@ 
        inits) /.
        { h[t0_] -> (solutionList[[maxOrder + 1]] /.
            locT -> t0),
          Derivative[n_][h][t0_] -> 
            (D[solutionList[[maxOrder + 1]], locT] /.
              locT -> t0)
        } // ExpandAll;
        
      (* Solve initials *)
      solutionList = solutionList /.
        Solve[Flatten[Map[Thread,
          Map[AsCoefficients[#, gaugeList, locEps]&, 
            subInits, {2}]]], constList][[1]] /. setGlobalRule;
       On[Solve::svars];
       solutionList
        
  ]

Finish

End[ ];

EndPackage[ ];