MultipleScales.mpl

# MODULE MultipleScales
#
# DESCRIPTION
# - A package used for solving systems of nonlinear partial differential
#   equations using the Method of Multiple Scales
#
# Written for Maple 2018
# Thomas Zdyrski (tzdyrski@ucsd.edu)
# 12 March 2018

MultipleScales := module()
description "Solves a system of nonlinear pdes using the Method of Multiple Scales":
option package:
local
  pdsolve,
  AppendSols,
  BothSides,
  ComplementarySols,
  DedupNamesFuncs,
  DedupRanking,
  DependentVars,
  DerivFuncs,
  Diff2Diffop,
  Diffop2Diff,
  ExtractCoeffs,
  EvalSolutions,
  HermitianInnerProduct,
  IndependentVars,
  MapRhs,
  RealEquations,
  RetainDerivs,
  SolveEquation,
  SolveEquations,
  SolveICs,
  SolveOrder,
  SolvingVars,
  ToAnalytic,
  TrigCoeff,
  TrigExpand,
  UniqueName,
  counter,
  UnapplyFuncs:

global
  `inttrans/expandc`:

export
  BasisCoeff,
  CombineSolutions,
  CompleteSystem,
  GenEqns,
  IndICs,
  IndFunc,
  IndFuncToName,
  IndName,
  IndGen,
  IndNameToFunc,
  SolveSystem,
  TruncateExpr,
  UnitTests,
  VerifyAndCombine,
  VerifySolution:

# Bug Fix submitted by vv at
# www.mapleprimes.com/questions/224385-Bug-Corrected-In-Inttranshilbert
`inttrans/expandc` := proc(expr, t)
local xpr, j, econst, op1, op2;
      xpr := expr;
      for j in indets(xpr,specfunc(`+`,exp)) do
          econst := select(type,op(j),('freeof')(t));
          if 0 < nops(econst) and econst <> 0 then
              xpr := subs(j = ('exp')(econst)*combine(j/('exp')(econst),exp),xpr)
          end if
      end do;
      for j in indets(xpr,{('cos')(linear(t)), ('sin')(linear(t))}) do
          if type(op(j),`+`) then
              op1:=select(has, op(j),t);
              op2:=op(j)-op1;
              if op(0,j) = sin then
                  xpr := subs(j = cos(op2)*sin(op1)+sin(op2)*cos(op1),xpr)
              else
                  xpr := subs(j = cos(op1)*cos(op2)-sin(op1)*sin(op2),xpr)
              end if
          end if
      end do;
      return xpr
end proc:

`pdsolve` := proc(eqn_sys)
local
  dep_vars,
  diff_ind_vars,
  doesnt_depend_on_diff_ind_vars,
  sols:
description "Overloads pdsolve to freeze all unknown functions that don't depend on the differentiated variables. This helps pdsolve solve equations it otherwise couldn't.":

  # Determine dependent variables to solve for
  if numelems({_rest}) >= 1 and type(_rest[1][1], {set,list}) then
    dep_vars := _rest[1][1]:
  else
    dep_vars := {}:
  end if:

  # Determine independent variables being differentiated
  convert(eqn_sys, Diff):
  indets(%, specop(`Diff`)):
  diff_ind_vars := map2(op,2,%):

  if diff_ind_vars = {} then
    # This isn't a differential equation; don't freeze variables

    sols := {solve(eqn_sys, dep_vars)}:

    # Check to make sure we found a solution
    if sols = {} then
      if _SolutionsMayBeLost=true then
        # Solve was just unable to solve system
        print("Warning: Unable to solve system; some solutions may be lost"):
        return NULL:
      else
        # System was inconsistent
        return NULL:
      end if:
    end if:

    return op(sols):
  end if:

  # Define type for unknown functions that don't depend on any of the
  # variables in diff_ind_vars
  doesnt_depend_on_diff_ind_vars :=
  And(
    unknown,
    function,
    satisfies(f -> {op(f)} intersect diff_ind_vars = {})
  ):

  # Check that system is consistent
  if numelems(eqn_sys) > 1 then
    frontend(PDEtools:-ConsistencyTest,[eqn_sys],
      [{`+`,`*`,`=`,`<>`,`<`,`<=`, Not(doesnt_depend_on_diff_ind_vars)}],
      dep_vars):
    if not % then:
      return NULL:
    end if:
  end if:

  # Only freeze those variables that don't depend on diff_ind_vars
  sols := {frontend(':-pdsolve', [eqn_sys],
    [{`+`,`*`,`=`,`<>`,`<`,`<=`, Not(doesnt_depend_on_diff_ind_vars)}],
    _rest)}:

  # Check to make sure we found a solution
  if sols = {} then
    # Solve was just unable to solve system
    print("Unable to solve system: solutions may be lost"):
    return NULL:
  end if:

  return op(sols):

end proc:

AppendSols := proc(
  existing_sols::Or(set,list)({name,function}=algebraic),
  new_sols::Or(set,list)({name,function}=algebraic),
  $)

description "Simply append new_sols after existing_sols":

  # Append `new_sols` to `existing_sols`
  [op(existing_sols), op(new_sols)]:

  return %:

end proc:

BasisCoeff := proc(
  expr::algebraic,
  basis_func::{function,constant},
  int_var::name:=NULL,
  {limits::algebraic:=infinity,
   regulate::truefalse:=true
  },
  $)
local
  new_basis_func,
  new_int_var,
  new_expr,
  integrand,
  regulator,
  normaliz:

description "Compute the coefficients of a basis function (or constant) `basis_func` in and expression `expr` using the inner product of the two functions. If `basis_func` depends on multiple variables, `int_var` must be the name of the variable over which to integrate. Integration limits are specified by `limits` (default: infinity) such that integration takes place from -`limits` to +`limits`. If regulate is true (default true), the integration will be regulated by a Gaussian regulator; may be slower.":

  # Define new variables that we can modify
  new_basis_func := basis_func:
  new_int_var := int_var:
  new_expr := expr:

  # Check if basis_func is a constant
  if type(basis_func, constant) and new_int_var = NULL then
    error "Must specify variable of integration `var_integration` if basis function is a constant.":
  end if:

  if new_int_var = NULL then
    # Determine integration variable as argument of basis_func
    new_int_var := op(basis_func);

    # It's possible the argument of basis_func is of the form c*f(x)
    # with c a constant; we only want f(x)
    indets(new_int_var, function):

    if % = {} then
    # It's possible the argument of basis_func is of the form c*x with
    # c a constant; we only want x
      indets(new_int_var, name):
    end if:

    # Convert set to series
    new_int_var := op(%):

    if numelems({%}) <> 1 then
      # Could not determine variable of integration
      error "Could not determine variable of integration; if either function depends on multiple variable, `int_var` must be passed":
    end if:
  end if:

  # Use HermitianInnerProduct to calculate dot product
  HermitianInnerProduct(Vector([expr]),Vector([basis_func]),new_int_var,':-limits'=limits):

  return %:

end proc:

BothSides := proc(
  Procedure::procedure,
  eqns::Or(set,list)({`=`,`<>`})
  )
description "Map `Procedure` to both sides of set\list `eqns`; pass any additioanl arguments to `Procedure`.":

  map2(map, Procedure, eqns, _rest):

  return %:

end proc:

ComplementarySols := proc(
  sols::{set,list}({name,function}=algebraic),
  existing_indets::{set,list}({name,function}),
  complementary_scales::{set,list}([{name,function},list(name)]),
  prev_sols::{set,list}({name,function}=algebraic),
  $)
local
  comp_name,
  comp_sol,
  new_indet,
  new_sols,
  new_existing_indets:

description "Given a set/list of solutions `sols` of the form {name,function}=algebraic, find any new indets not present in the set/list `existing_indets`. Then, for each name/function in `complementary_scales`, check if any of the new indets are present in that term's solution. If so, replace the indet with a new function with restricted operators as supplied by the second element of the complementary_scales list. The set/list of previous solutions `prev_sols` is used to determine which new indets are complementary solutions and which are not.":

  new_sols:= sols:

  # "Complementary solutions" are a bit tricky to define for PDEs. We
  # will only consider new variables appearing solely in solutions for
  # new variables as possible complementary solutions. Eg if a solution
  # contains a term like `_D32(t__1,t__2) = _D33(t__2)`, we *won't*
  # consider `_D33` at complementary solution because `_D32` appeared in
  # previous orders' solutions. Another eg: if at order `n` the solution
  # contains `y_n(t) = _D1*sin(t)`, we will consider
  # `_D1` a complementary solution (if it doesn't appear in the solution
  # of any other terms that appeared in previous solutions) because
  # `y_n(t)` has not appeared in any previous solutions. One more eg: if
  # previous solutions were `_D12(t__0,t__1,t__2) = _D13(t__1,t__2)` and
  # current solution is `y__3(t)=_D14*sin(t), phi_3(t)=_D15*cos(t),
  # _D13(t__1,t__2)=_D16(t__2)`, then `_D14` and `_D15` would be
  # complementary solutions since neither `y__3(t)` nor `phi__3(t)`
  # appeared in previous solutions, but `_D16(t__2)` would not be a
  # compementary solution since `_D13` appeared in a previous solution.
  # This criterion attempts to capture the intuitive concept of
  # "complementary solutions" formally.

  # Add terms appearing as solutions to previous solutions to set of
  # existing_indets
  eval(prev_sols, sols):
  indets(%, assignable):
  new_existing_indets := % union convert(existing_indets,set):

  # for each complementary name
  unassign('comp_name'):
  for comp_sol in complementary_scales do

    # Get solution for comp_name
    eval(op(1,comp_sol), new_sols):

    # Get new indets
    indets(%, 'assignable') minus new_existing_indets:

    for new_indet in % do
      # Define new operands
      select(member,[op(new_indet)], op(2,comp_sol)):

      # If no operands, make constant
      if numelems(%) = 0 then
        new_sols := eval(new_sols, op(0,new_indet) =
          unapply(UniqueName(),op(new_indet))):
      else
        new_sols := eval(new_sols, op(0,new_indet) =
          unapply(UniqueName()(op(%)),op(new_indet))):
      end if:

    end do:

  end do:

  return new_sols:

end proc:

DedupNamesFuncs := proc(
  funcs::Or(set,list)({name,function}),
  {ind_vars::Or(set,list)({name,function}):={}},
  $)
description "Given a set/list of functions, remove duplicates; optionally, remove functions with arguments not included in the set/list ind_vars.":

  # Convert to list
  convert(funcs, list):

  if numelems(ind_vars) > 0 then
    # Remove functions depending on arbitrary independent variables
    remove((x,set) -> type(x,function) and not({op(x)} subset set ), %, convert(ind_vars,set)):
  end if:

  # Remove duplicate functions
  ListTools:-MakeUnique(%, 1,
    proc(x,y) type(x,'function') and type(y, 'function') and
      evalb(op(0,x) = op(0,y)) end proc):

  # Remove duplicate names
  ListTools:-MakeUnique(%, 1,
    proc(x,y) type(x,'name') and type(y, 'name') and
      evalb(x = y) end proc):

  # Convert back to original type
  convert(%, whattype(funcs)):

end proc:

DedupRanking := proc(
  ranking::list(Or(set,list)),
  $)
local
  rev_ranking,
  new_ranking,
  j:
description "Remove duplicates from lower order ranking.":

  # Reverse ranking
  rev_ranking := ListTools:-Reverse(ranking):

  new_ranking := []:
  unassign('j'):
  # Remove duplicates from lower order in the rev_ranking
  for j from 1 to numelems(rev_ranking) do
    # Remove indets from current element that appear in previous
    # elements
    remove(member, rev_ranking[j], `union`(op([1..j-1], map(convert,
      rev_ranking, set)))):

    if % <> {} and % <> [] then
      # Append to new_rev_ranking
      new_ranking := [op(new_ranking), %]:
    end if:
  end do:

  # Undo ranking reversal
  ListTools:-Reverse(new_ranking):

  return %:

end proc:

DependentVars := proc(
  exprs::Or(set,list),
  {ind_vars::Or(set,list)({name,function}):={}},
  $)
local
  arguments:

description "Return a set of variables which appear as functions in set of equations or inequations `eqn_sys`. If `ind_vars` is passed, only functions depending on a subset of `ind_var` will be returned.":

  # Select all functions
  indets(convert(exprs,Diff),And(assignable,function,unknown)):

  # Convert to list
  convert(%, list):

  # Remove duplicate functions
  DedupNamesFuncs(%, ':-ind_vars'=ind_vars):

  # Convert to set
  convert(%, set):

  return %:

end proc:

EvalSolutions := proc(
  exprs,
  solutions::Or(set,list)({name,function}=algebraic),
  only_derivs::Or(set,list)({name,function}):={},
  $)
local
  eval_fcns,
  exclude_fcns:

description "Evaluates `exprs` using `solutions`; any functions in `only_derivatives` will only have their derivatives evaluated, but their value will not be substituted otherwise":

  # Remove any `only_derivs` from `solutions`
  exclude_fcns, eval_fcns := selectremove((expr,arg) ->
    member(lhs(expr),arg), convert(solutions,D), only_derivs):

  # Unapply solutions to get them in functional form
  eval_fcns := UnapplyFuncs(eval_fcns):
  exclude_fcns := UnapplyFuncs(exclude_fcns):

  # Insert `solutions` to be evaluated
  eval['recurse'](convert(exprs,D), eval_fcns):

  # Evaluate derivatives
  value(%):

  # Evaluate derivatives of `only_derivs` functions
  evalindets(convert(%,D),
    typefunc(
      typefunc(
        {map(identical,map2(op,0,convert(only_derivs,set)))},
        Or(
          identical(D),
          specindex(posint,D),
          And(
            typefunc(
              identical(`@@`)
            ),
            anyfunc(identical(D), posint)
          )
        )
      )
    ), eval['recurse'], [op(exclude_fcns), op(eval_fcns)]):

  return %:

end proc:

DerivFuncs := proc(
  eqns::Or(set,list)(equation),
  deriv_var::name,
  $)
local
  deriv_wrt_var:

description "Given a set/list of equations `eqns`, return a set/list of sets of all functions that appear in `eqns` and are differentiated (at least once) with respect to the name `deriv_var`.":

  # Define depends_on_deriv_var
  deriv_wrt_var := And(
    typefunc(identical(diff)),
    anyfunc(function,identical(deriv_var))
  ):

  # Expand exprs to ensure derivatives appear as monomials
  expand(eqns):

  # Convert to diff-notation
  convert(%, diff):

  # Get all derivatives wrt deriv_var
  map(indets, %, deriv_wrt_var):

  # Extract functions from derivatives
  map(eq -> indets(`union`(eq),assignable(function)), %):

  return %:

end proc:

Diff2Diffop := proc(
  eqns::{set,list}({equation,algebraic}),
  {vars::{set,list}:={},
   funcs::{set,list}:={}
   },
  $)
local
  ind_vars,
  repl_funcs,
  replaced_eqns,
  prev_replaced_eqns,
  diff_vars,
  var:

description "Given a set/list of equations or expressions `eqns`, returns the equations with all derivative replaced by symbols Dname with name representing the differentiating variable. Also returns a set of equations of the form Dname=name for each differentiation variable. By default, all independent variables are used; if only some should be transformed, pass a set/list `vars` of differentiation variable names to transform. If a set/list of functions `func` is passed (eg {f(x),g(t,x)}), only derivatves of those functions will be transformed.":

  if vars <> {} then
    ind_vars := vars:
  else
    # Get list of variables
    ind_vars := IndependentVars(eqns):
  end if:

  # Specify which functions to replace
  if numelems(funcs) > 0 then
    repl_funcs := Or(op(map(identical,funcs))):
  else
    repl_funcs := function:
  end if:

  # Create empty variables to store results
  replaced_eqns := expand(convert(eqns,diff)):
  prev_replaced_eqns := replaced_eqns:

  diff_vars := {}:

  # For each independent variable, replace
  for var in ind_vars do

    # Loop until applying the transform no longer changes anything (this
    # is necessary to handle multiple derivatives
    do
      # Replace derivative
      replaced_eqns := evalindets(replaced_eqns,
        And(
          anyop(repl_funcs,identical(var)),
          typefunc(identical('diff'))
        ),
        f -> cat(`D`,var)*op(1,f)):

      # Check if the result has stopped changing
      if prev_replaced_eqns = replaced_eqns then
        break:
      else
        prev_replaced_eqns := replaced_eqns:
      end if:
    end do:

    # Add differential variable to diff_vars
    diff_vars := diff_vars union {cat(`D`,var) = var}:

  end do:

  return replaced_eqns, diff_vars:

end proc:

Diffop2Diff := proc(
  eqns::{set,list}({equation,algebraic}),
  diff_vars::{set,list}(name=name),
  $)
local
  replaced_eqns,
  prev_replaced_eqns,
  diff_repl:

description "Given a set/list of equations or expressions `eqns` and a set/list of differentiation variables `diff_vars` of the form Dname=name, replaces each instance of Dname with diff( ,name).":

  # Note: it is important here that, when converting back to
  # derivatives, we take the derivative of *all* terms.
  # For instance, exp(t)*diff(f(t),t) converts to Dt*exp(t)*f(t) and converts
  # back to diff(exp(t)*f(t),t); this may seem wrong at first, but
  # Diffop2Diff is used when taking transposes. the transpose of
  # exp(t)*diff(<something>,x) is in fact -diff(exp(t)*<something>,x),
  # as can be seen by integrating by parts.

  # Create empty variables to store results
  replaced_eqns := expand(eqns):
  prev_replaced_eqns := replaced_eqns:

  # For each independent variable, replace
  for diff_repl in diff_vars do
    # Loop until applying the transform no longer changes anything (this
    # is necessary since powers (eg Dname^3) will only trigger this
    # condition once
    do
      # Replace derivative operators
      replaced_eqns := evalindets(replaced_eqns,
        And(
          satisfies(f -> has(f,lhs(diff_repl))),
          specop(`*`)
        ),
        f -> diff(f/lhs(diff_repl),rhs(diff_repl)));
      # Check if the result has stopped changing
      if prev_replaced_eqns = replaced_eqns then
        break:
      else
        prev_replaced_eqns := replaced_eqns:
      end if:
    end do:
  end do:

  return replaced_eqns:

end proc:

ExtractCoeffs := proc (
  exprs::Or(set,list)(equation),
  orth_basis::Or(set,list)(function),
  {exclude_vars::Or(set,list):={},
   extract_type::{"trig","innerProduct"}:="trig",
   limits::algebraic:=Pi},
  $)

description "Given a set or list of functions `orth_basis`, extract the components of each basis function from each element of set/list `exprs`.  If `extract_type` is 'trig' (default), then `orth_basis` should be trigonometric functions and `exclude_vars` should consist of the arguments to `orth_basis`. If `extract_type` is `innerProduct`, a hermitian inner product is used to calculate the coefficients of `orth_basis` and `limits` (default Pi) are the +/- limits of integration; note: `innerProduct` is much slower than `trig`.":

  # Isolate coefficients of basis functions
  if extract_type = "trig" then
    map[3](op@BothSides, TrigCoeff, exprs, orth_basis, exclude_vars):
  else
    map[3](op@BothSides, BasisCoeff, exprs, orth_basis, ':-limits'=limits):
  end if:

  return %:

end proc:

HermitianInnerProduct := proc(
  bra::Vector,
  ket::Vector,
  int_var::name,
  {limits::algebraic:=infinity,
   op_list::Or(set,list)(name=name):={}},
  $)
local
  additional_op_list,
  identical_diffop,
  contains_diffop,
  assumptions,
  new_bra:

description "Given a two Vectors `bra` (possibly containing differential operators) and `ket` (not containing differential operators), calculate the hermitian inner product by integrating over variable `int_var`. If algebraic `limits` is provided, then the integral has limits of -limits..limits. If a set/list `op_list` of differentiation variables `diff_vars` of the form Dname=name, replaces each instance of Dname with diff( ,name).":

  # Get list of all dependent and independent variables (so we can
  # assume real) except derivatives
  `union`(indets(convert(bra,list),assignable),
    indets(convert(ket,list),assignable)):
  % minus convert(map(lhs,op_list),set):

  # Assume real
  map(x -> x::real, %):

  # Assume derivatives are anti-Hermitian
  assumptions := % union map(x -> x::imaginary, map(lhs,op_list)):

  ## We want to ensure all derivatives are conjugated correctly.
  # Convert derivatives to additional_op_list
  new_bra, additional_op_list := Diff2Diffop(convert(bra,list)):

  # Add new additional_op_list to assumptions
  assumptions union map(x -> x::imaginary, map(lhs,additional_op_list)):

  # Take conjugate of bra to convert vector to 1-form
  conjugate(new_bra) assuming op(%):

  # Convert back to derivatives
  new_bra := Vector(Diffop2Diff(%, additional_op_list)):

  # Element-wise multiply each element of conjugate(ket) and bra
  # without taking the conjugate (we already did that)
  LinearAlgebra:-DotProduct(new_bra,ket, conjugate=false):

  # Check if there are diff_ops
  if numelems(op_list) <> 0 then
    # Generate type containing diffops
    identical_diffop := Or(op((map(identical, map(lhs,op_list))))):
    contains_diffop := satisfies(x -> hastype(x, identical_diffop)):

    # Freeze all non-diffop terms
    expand(new_bra):
    new_bra := map2(map, evalindets, %, And(Not(contains_diffop),
      function), x->freeze(x)):

    # Element-wise multiply each element of conjugate(ket) and bra
    # without taking the conjugate (we already did that)
    LinearAlgebra:-DotProduct(new_bra,ket, conjugate=false):

    # Convert back to diff-form
    Diffop2Diff({%}, op_list)[1]:

    # Expand
    expand(%):

    # Thaw
    thaw(%):
    do
      # Loop because it's possible frozen terms contain other frozen terms
      thaw(%):
      if % = %% then
        break:
      end if:
    end do:
  end if:

  # Convert trig to exp to help integration
  evalindets(%, And(function,satisfies(f -> has(f,int_var))), convert,
    exp) assuming op(assumptions);
  expand(%):

  # Replace and Dirac with Dirac*lim so that they don't vanish when
  # taking int / lim
  eval(%, Dirac= (proc()  lim*Dirac(args) end proc));

  # Calculate inner product
  int(%,int_var=-lim..lim):

  # Take integration limits to `limits` and divide by `limits`
  limit(%/lim,lim=limits) assuming real:

  if has(%, int) then
    print("Warning: Integral in HermitianInnerProduct couldn't fully evaluate; some solutions may be lost."):
    return false:
  end if:

  return %:

end proc:

IndependentVars := proc(
  exprs::Or(set,list),
  $)
local
  arguments:

description "Return a set of variables which appear as parameters of functions in set of equations or inequations `eqn_sys`":

  # Select all functions
  indets(exprs,function):

  # It's possible we have functions of functions of ... of arguments;
  # recursively select arguments the list of arguments doesn't change
  arguments := %:
  while 1 = 1 do
    map(op, arguments):

    if % = arguments then
      arguments := %:
      break:
    end if:

    arguments := %:
  end do:


  # It's possible we have some constants now; remove them
  indets(%, assignable):

  return %:

end proc:

MapRhs := proc (
  fcn::Or(procedure,name),
  equn::Or(equation,{set,list}(equation))
  )
local
  rest_args:
description "Maps function `fcn` to the right-hand-side of `equn`":

  # It seems the special variable `_rest` can't be accessed inside a
  # `map` command; instead, store the result to a temporary variable
  # `rest_args` that can be accessed with `map`
  rest_args := _rest:

  map(expr -> lhs(expr) = fcn(rhs(expr), rest_args), equn):

end proc:

SolveOrder := proc(
  PDEs::{set,list}(Or(`=`, `<>`)),
  solving_order::integer,
  mult_scales::list(name),
  {const::Or(set,list)(name):={},
   prev_sols::Or(set,list)({name,function}=algebraic):={},
   exclude_vars::Or(set,list)({name,function}):={},
   ind_vars::Or(set,list)({name,function}):={},
   constraints::Or(set,list)({`=`({name,function},constant),`<>`({name,function},constant),
     `<`({name,function},constant),`>`({name,function},constant),
     `<=`({name,function},constant),`>=`({name,function},constant)}):={},
   extract_type::{"trig","innerProduct"}:="trig",
   secular_allowed::Or(set,list)(function):={},
   complementary_scales::{set,list}([{name,function},list(name)]):={}},
  $)
local
  eval_pdes,
  diff_funcs,
  op_pdes,
  op_list,
  lin_pdes,
  nonlin_pdes,
  hom_sols,
  new_unknowns,
  compatibility,
  term,
  zero_eigenvector,
  compat_sols,
  solving_system,
  existing_indets,
  hom_terms,
  inhom_terms,
  unknowns,
  hom_eqns,
  null_space_vectors,
  new_ind_vars,
  HomInhomTerms,
  ConjTransHomEqns,
  NullSpaceVecs,
  GenAndSolveSys:

description "Solve a set\list of partial differential equations `PDEs` using a set\list of previous solutions `prev_sols` of the form `function(args) = expression`. `solving_order` denotes which order in the perturbative parameter `PDEs` satisfies. `mult_scales` is a list of names of independent variables; secular terms in these variables will be avoided (in order) unless the term becoming secular appears in `secular_allowed`. `const` is a set of constants that should not be solved for (arbitrary constants of the problem).  `exclude_vars` are functions/variables whose values should not be substituted into the answer (though derivative of these quantities will still be evaluated).  Additionally, `constraints` are any additional equation/inequalities of the form 'function/symbol =/</>/<=/>=/<> constant'. `extract_type` is passed to ExtractCoeffs; see ExtractCoeffs for a description. If `ind_vars` is supplied, only functions of these variables will be considered. If a set/list of lists `complementary_scales` of the form [{name,function},[names]] is supplied, then any function appearing in the first element (either a name or function) will only have homogenous terms depending on the variables `names`. This is useful for restricting higher-order terms from having too many free parameters.":

HomInhomTerms := proc(
  eqns::Or(set,list)(equation),
  funcs::Or(set,list)(function),
  $)
local
  unique_const,
  contains_funcs,
  hom_terms,
  inhom_terms,
  op_eqns,
  common_factor:

description "Given a set/list of differential equations `eqns` and a set/list of functions `funcs`, return a set/list of the terms depending on `funcs` followed by a set/list of the terms that don't.":

  # Hack-y: generate a random constant to add to each equation so we can
  # guarantee every equation is a sum and not a monomial; add two
  # constants together in case the equation was 0=0, we need to to still
  # get a sum
  UniqueName():
  unique_const := % + UniqueName():

  # Define proc that returns true for terms linear in funcs after
  # converting all diffs to diff operators; also return true for
  # independent variables
  contains_funcs :=f -> hastype(
      f, Or(op(map(identical, funcs)))):

  # Move all terms to the same side
  map(x->lhs(x)-rhs(x), eqns):

  # Make sure we have a sum of terms
  expand(%):
  map(`+`, %, unique_const):

  # Convert to diff
  convert(%, diff):

  # Separate terms linear in funcs
  map(f -> [selectremove(contains_funcs, f)], %):
  hom_terms, inhom_terms := map2(op,1,%), map2(op,2,%):
  inhom_terms := map(`-`, inhom_terms, unique_const):

  return hom_terms, inhom_terms:

end proc:

ConjTransHomEqns := proc(
  hom_eqns::list(algebraic = 0),
  funcs::Or(set,list)(function),
  $)
local
  op_eqns,
  op_list,
  diffop_names,
  linear_op,
  nonlinear_op,
  null_space_ind_vars,
  null_space_vector,
  hom_eq,
  identical_diffop,
  contains_diffop,
  contains_null_vector:

description "Must pass a list of homogenous equations (algebraic = 0) that are linear in all the functions in the set/list funcs. Returns a list of functions where the linear operator has been replaced by its conjugate transpose and applied to a new set of functions; the second returned parameter is the list of these new functions.":

  # Convert equations to diff form
  op_eqns, op_list := Diff2Diffop(hom_eqns):

  # Get list of diffops
  diffop_names := map(lhs,op_list):

  # Generate type containing diffops
  identical_diffop := Or(op((map(identical, diffop_names)))):
  contains_diffop := satisfies(x -> hastype(x, identical_diffop)):

  # Convert to a matrix equation
  linear_op, nonlinear_op := LinearAlgebra:-GenerateMatrix(op_eqns, funcs):

  if not LinearAlgebra:-Equal(nonlinear_op,
    Vector(LinearAlgebra:-Dimension(nonlinear_op),fill=0)) then
      error(cat("Nonlinearity in order-",solving_order,", ", mult_scales[1],
      "-derivative equation: ", op_eqns, "\nThis is not yet implemented.")):
  end if:

  # Generate test functions to determine null space of Dagger(lin_pdes)
  null_space_ind_vars := IndependentVars(hom_eqns) intersect ind_vars:

  unassign('j'):
  null_space_vector := [seq(UniqueName()(op(null_space_ind_vars)),
    j=1..LinearAlgebra:-RowDimension(linear_op))]:

  # Create type
  contains_null_vector := Or(op(map(identical,null_space_vector))):

  # Assume all quantites (except derivatives) are real
  indets(linear_op,assignable) union indets(null_space_vector,assignable):
  % minus diffop_names:
  map(x -> x::real, %):
  # Assume derivatives are anti-Hermitian
  % union map(x -> x::imaginary, diffop_names):

  # Regenerate linear, homogenous part assuming real
  Physics:-Dagger(linear_op).Vector(null_space_vector) assuming op(%):
  convert(%, list):

  # Convert back to diff-form
  Diffop2Diff(%, op_list):

  # Set equal to zero
  map(x -> x=0, %):

  return %, null_space_vector:

end proc:

NullSpaceVecs := proc(
  hom_eqns::Or(set,list)(algebraic = 0),
  funcs::Or(set,list)(function),
  diff_var::name,
  {const::Or(set,list)(name):={},
   ind_vars::Or(set,list)({name,function}):={}},
  $)
local
  trans_eqns,
  trans_funcs,
  trans_sols,
  trans_sol,
  null_vectors,
  all_null_vectors,
  op_eqns,
  op_list,
  op_sols,
  dname,
  diffop,
  max_degree,
  ranking,
  j,
  replacements,
  lim:

description "Given a set/list of homogenous differential equations and a set/list of functions, returns a set of lists representing a basis for the null space of the conjugate of the differntial equation's operator.  A name `diffvar` should be passed to indicate the variable wrt to which derivatives are being considered. Sets/lists of names (`const`) and names/functions (`ind_vars`) are passed to SolveEquations when calculating the null space. If any null vectors depend on diff_var, they will be replaced by a Dirac delta function Dirac(diff_var - diff_var__dummy).":

  # Check that our homogenous equation system isn't just an array of
  # 0=0.
  if convert(hom_eqns,set) = {0=0} then
    return {}:
  end if:

  # Take conjugate transpose of linear operator
  trans_eqns, trans_funcs := ConjTransHomEqns(hom_eqns, funcs):

  # Rank solving_vars so that terms with more derivatives are solve
  # first; we want solutions like _D1(t) = diff(_D2(t),t), not _D2(t) =
  # int(_D1(t),t)
  op_eqns, op_list := Diff2Diffop(trans_eqns, ':-funcs'=trans_funcs):

  # We will count all derivatives equally, so replace them all with a
  # temporary name
  dname := UniqueName():
  for diffop in map(lhs, op_list) do
    op_eqns := eval(op_eqns, diffop = dname):
  end do:

  # Move all terms to one side
  op_eqns := map(x -> lhs(x) - rhs(x), op_eqns):

  # Determine highest power of dname
  map(degree, op_eqns, dname):
  max_degree := max(%):

  unassign('j'):
  ranking := []:

  # Put higher degrees earlier in the ranking (ie later in the list)
  for j from max_degree by -1 to 0 do
    # Collect all terms of order-j in dname
    map(coeff, op_eqns, dname, j):

    # Convert to set
    convert(%, set):

    # We only care about terms, not constant coefficients; remove them
    indets(%, assignable) minus convert(const,set) minus
      convert(ind_vars,set):

    # Only consider functions that are being differentiated
    % intersect convert(trans_funcs,set):

    # Add to ranking
    ranking := [%, op(ranking)]:
  end do:

  # Dedup ranking (in reverse since we want the *higher* ranking
  # duplicates removed
  ListTools:-Reverse(ranking):
  DedupRanking(%):
  ranking := ListTools:-Reverse(%):

  # Get list of unknowns
  unknowns := indets(trans_eqns, 'assignable'):

  # Solve homogenous case
  trans_sols := SolveEquations(trans_eqns, ':-const'=const,
    ':-ranking'=ranking, ':-ind_vars'=ind_vars):

  # Only keep solutions that are fully solved
  trans_sols := remove(has, %, int):

  # Check that we generated null vectors
  if % = {} then
    print("Warning: Couldn't compute null space; some solutions may be lost."):
    return {}, op_list:
  end if:

  # Get null vectors for each solution set
  all_null_vectors := {}:
  unassign('trans_sol'):
  for trans_sol in trans_sols do

    # Determine new undetermined coefficients
    map(rhs,trans_sol):
    new_unknowns := indets(%, 'assignable'):

    # Remove previously known quanitites
    new_unknowns :=  new_unknowns minus unknowns minus map(lhs,op_list):

    # Also add `free` parameters of the form  X=X that can take any value
    select(x -> evalb(lhs(x) = rhs(x)), trans_sol):
    map(rhs, %):
    new_unknowns := new_unknowns union %:

    # Also add functions that trans_eqns didn't depend on (so they're
    # allowed to be anything)
    new_unknowns := new_unknowns union indets(trans_funcs, 'assignable')
      minus convert(ind_vars,set):

    # Generate null vectors
    null_vectors := {}:
    unassign('term'):
    for term in new_unknowns do
      # Set all other coefficients equal to zero
      map(x -> x=lim, new_unknowns minus {term}):

      if depends(term, diff_var) then
        # Set current coefficient to Dirac delta function
        replacements := % union {term = Dirac(diff_var -
          cat(diff_var,`__dummy`))}:
      else
        # Set current coefficent to unity
        replacements := % union {term = 1}:
      end if:

      # Evaluate solutions with only current coefficient non-zero
      map(x -> lhs(x) = eval(rhs(x), replacements), trans_sol):

      # Put the solutions in the same ordering as the trans_funcs vector
      eval(trans_funcs, %):

      # In case some functions weren't included in tran_equations (and
      # hence are free), also replace tran_funcs with replacements
      eval(%, replacements):

      # Take lim to 0
      map(limit, %, lim=0,right):

      # Add to the set of null_vectors
      null_vectors := null_vectors union {%}:

    end do:

    # Add null_vectors to set of all null vectors
    all_null_vectors := all_null_vectors union null_vectors:

  end do:

  return all_null_vectors, op_list:

end proc:





  # Insert previous solutions
  EvalSolutions(PDEs, prev_sols, exclude_vars):

  # Convert to list (order will matter)
  eval_pdes := convert(%, list):

  # Get a list of functions differentiated wrt mult_scales[1]
  DerivFuncs(%, mult_scales[1]):

  # Combine functions from each equation
  `union`(op(%)):

  # Remove `allowed_secular` terms
  diff_funcs := % minus secular_allowed:

  # Convert to list because ordering will be important
  diff_funcs := convert(diff_funcs, list):

  # Check that we don't have more functions than equations
  if numelems(%) > numelems(eval_pdes) then
    error cat("Warning: underdetermined system at order ", solving_order):
  end if:

  # Define a new set of ind_vars
  new_ind_vars := ind_vars:

if diff_funcs <> [] then
  # Split into homogenous and inhomgenous parts
  hom_terms, inhom_terms := HomInhomTerms(eval_pdes, diff_funcs):
else
  diff_funcs := indets(eval_pdes, 'assignable') minus const minus
    ind_vars:
  inhom_terms := {0}:
end if:

# See if there are any inhomogenous terms
if convert(inhom_terms,set) <> {0} then
  # There ar inhomogenous terms; generate compatibility conditions

  # Generate equations from linear terms
  convert(hom_terms, list):
  hom_eqns := map(x -> x=0, %):

  # Generate list of null vectors
  null_space_vectors, op_list := NullSpaceVecs(hom_eqns, diff_funcs,
    mult_scales[1], ':-const'=const, ':-ind_vars'=new_ind_vars):

  # Check that we generated null vectors
  if null_space_vectors = {} then return {} end if:

  # Add dummy variable (possibly created by NullSpaceVecs) to set of
  # ind_vars
  new_ind_vars := {op(ind_vars), cat(mult_scales[1],`__dummy`)}:

  # Create compatibility conditions
  compatibility := {}:
  for zero_eigenvector in null_space_vectors do

    # Form product <0-eigenvector | non-linear equations>
    HermitianInnerProduct(Vector(zero_eigenvector),Vector(inhom_terms),
      mult_scales[1], ':-limits'=infinity, ':-op_list'=op_list):

    # If HermitianInnerProduct couldn't evaluate integral, return an
    # empty set to drop this solution
    if % = false then return {} end if:

    # Add to compatibility set
    compatibility := compatibility union {%=0}:

  end do:

  # Remove 0=0 equation
  compatibility := compatibility minus {0=0}:

  compat_sols := {}:
  if numelems(mult_scales) > 1 and compatibility <> {} then
  # There are more scales and equations to consider

    # Apply same algorithm to compatibilty conditions
      compat_sols := SolveOrder(
        compatibility, solving_order, mult_scales[2..],
        ':-const'=const,
         ':-prev_sols'=prev_sols,
         ':-exclude_vars'=exclude_vars,
         ':-ind_vars'=new_ind_vars,
         ':-constraints'=constraints,
         ':-extract_type'=extract_type,
         ':-secular_allowed'=secular_allowed,
         ':-complementary_scales'=complementary_scales):

    # If no solution found, return empty solution set
    if % = {} then return {} end if:

  elif compatibility <> {} then
  # There are more equations to consider (but *not* more
  # scales)
    compat_sols := SolveEquations(compatibility, ':-const'=const,
      ':-ind_vars'=new_ind_vars):

    # If no solution found, return empty solution set
    if % = {} then return {} end if:

  end if:
else
  # There aren't any compatibility conditions
  compat_sols := {}:
end if:

  # Add constraints equations to system of pdes
  solving_system := [op(eval_pdes), op(EvalSolutions(constraints,
    prev_sols, exclude_vars))]:

  # Determine existing indets
  existing_indets := indets(solving_system, 'assignable'):

  if compat_sols <> {} then
    # Define proc to make steps clearer: Note, this proc cannot be moved
    # outside of this if/then statement since it relies on externally
    # defined variables solving_system, exclude_vars, const, ind_vars,
    # and diff_funcs
    GenAndSolveSys := proc(
      sol::set({name,function}=algebraic),
      $)
    local
      solving_vars:

      # Add solutions of compatibility equations to `solving_vars`
      solving_vars :=convert(diff_funcs,set) union map(lhs,sol):

      # Evaluate equations using compatibility solutions; each set of
      # compatibility solutions should be evaluated separately
      EvalSolutions(solving_system,sol,exclude_vars):

      # Solve
      SolveEquations(%, ':-const'=const, ':-ind_vars'=new_ind_vars,
          ':-ranking'=[solving_vars]):

      # Add compatibility solutions
      map(`union`, %, sol):

      return %:
    end proc:

    # Solve the set of pdes corresponding to each possible
    # compatibility-equation solution set
    map(GenAndSolveSys, compat_sols):
  else

    # Just solve `solving_system`
    {SolveEquations(solving_system, ':-const'=const,
      ':-ranking'=[convert(diff_funcs,set)],
      ':-ind_vars'=new_ind_vars)}:

  end if:

  # Only keep solutions that are fully solved
  remove(has, %, int):

  # Package all solutions to various compatibility equations together
  map(op,%):

  # Restrict complementary solutions
  map(ComplementarySols, %, existing_indets, complementary_scales, prev_sols):

  return %:

end proc:

RealEquations := proc(
  eqns::Or(set,list)(equation),
  siderels::Or(set,list)({name,function}=algebraic):={},
  {assumptions::Or(set,list)(Or(name::property,And(relation,anyop(name,algebraic)))):={}},
  $)
local
  ind_vars:

description "Given a set/list of equations, it verifies that all equations are real when their names/functions are assumed positive. If `siderels` are supplied, they are used to evaluate `eqns` before checking realness. If `assumptions` are supplied, they will be applied while evaluating eqns.":

  # Determine independent variables
  ind_vars := IndependentVars(eqns):

  # Evaluate eqns using siderels
  EvalSolutions(eqns, siderels):

  # Evaluate lhs and rhs of each equation assuming all functions/names
  # are positive
  map(x -> [evalc(Im(lhs(x))), evalc(Im(rhs(x)))], %) assuming op(assumptions), positive:

  # Check that each side was zero
  map(x->verify(x,[0,0]), %):

  # Check that all equations were real
  `and`(op(%)):

  return %:

end proc:

RetainDerivs := proc(
  expr::algebraic,
  deriv_var::name,
  $)
local
  deriv_wrt_var,
  depends_on_diff:

description "Given an expression `expr`, remove all terms that do not contain a derivative wrt `deriv_var`.":

  # Define depends_on_deriv_var
#  depends_on_deriv_var := satisfies(f -> `or`(op(map(type,{op(f)},identical(deriv_var))))):
#  depends_on_deriv_var := anyfunc(function,identical(deriv_var)):
#  deriv_wrt_var := And(typefunc(identical(diff)), depends_on_deriv_var):
  deriv_wrt_var := And(
    typefunc(identical(diff)),
    anyfunc(function,identical(deriv_var))
  ):
  depends_on_diff := satisfies(f -> hastype(f, deriv_wrt_var)):

  # Expand exprs to ensure derivatives appear as monomials
#  Physics:-Expand(expr):
  expand(expr):

  # Convert to diff-notation
  convert(%, diff):

  # Check if we have a sum of terms
  if type(%,`+`) then
    # Select terms that depend on diff wrt deriv_var
    select(type, %, depends_on_diff):
    if %=0 then
      # No terms found
      return:
    else
      return %:
    end if:
  else
    # Only a single monomial; check entire thing
    if type(%, depends_on_diff) then
      return expr:
    else
      # Doesn't depend on diff wrt deriv_var; return nothing
      return:
    end if:
  end if:

end proc:

SolveEquations := proc(
  eqn_sys::Or(set,list)({`=`,`<>`}),
  {const::Or(set,list)(name):={},
   ranking::list(Or(set,list)({name,function})):=[],
   ind_vars::Or(set,list)({name,function}):={}},
  $)
local
  split_sys,
  eqns,
  results:
description "Split system and apply SolveEquation to each system.":

  # Split system
  convert(eqn_sys,set) minus {0=0}:
  expand(%):

  split_sys := {PDEtools:-splitsys(%,
    indets(convert(eqn_sys,diff), assignable) minus convert(const,set)
    minus convert(ind_vars,set))}:

  # Remove equations of the form f(x) = f(x)
  split_sys := remove(type, split_sys, And(set(equation), satisfies(x ->
    numelems(x) = 1), satisfies(x -> lhs(op(x)) = rhs(op(x))))):

  # If the constraints decouple completely, remove them as well
  split_sys := remove(type, split_sys, Or(`<>`,
    {set,list}(`<>`))):

  results := {{}}:
  unassign('eqns'):
  for eqns in split_sys do
    SolveEquation(eqns, args[2..]):

    # Remove solutions of the form f(x) = f(x)
    map2(remove, type, %, satisfies(f -> lhs(f) = rhs(f))):

    # Append to previous results
    results := map(x -> op(map(`union`, results, x)), %):

  end do:

  return results:

end proc:

SolveEquation := proc(
  eqn_sys::Or(set,list)({`=`,`<>`}),
  {const::Or(set,list)(name):={},
   ranking::list(Or(set,list)({name,function})):=[],
   ind_vars::Or(set,list)({name,function}):={}},
  $)

local
  new_eqn_sys,
  new_ranking,
  existing_names,
  expr,
  new_solving_vars,
  orig_vars,
  params,
  pde_equations,
  constraints,
  solving_system,
  other_funcs,
  X,
  _A,
  _phi,
  new_eq,
  denom_eq,
  sols,
  new_sols,
  sol,
  temp_sol,
  unique_const,
  repl_name,
  CheckConstraints:

description "Solve a set of equation/inequations by determining the depenendent variables (excluding those in `consts`) and using Solve.  `ranking` is a list of sets\lists that gives the solving_vars of solving variables. If `ind_vars` is supplied, only functions of these variables will be considered.":

CheckConstraints := proc(
  sols::Or(set,list)({name,function}=algebraic),
  constraints::Or(set,list)({equation,`<>`}),
  $)
description "Given a set/list of solutions, check if they can satisfy a given set/list of constraints. If they can, return the solutions; otherwise, return NULL.":

  # Apply to constraints
  eval(constraints, sols):

  # See if constraints could be satisfied
  `and`(op(map(coulditbe, %))):

  # If true, return sols; else, return NULL
  if % then
    return sols:
  else
    return:
  end if:
end proc:

  new_eqn_sys := {}:
  sols := {}:

  # Determine all functions/names
  existing_names := indets({op(eqn_sys), op(ranking), op(const),
    op(ind_vars)}, 'assignable'):

  # Separate equations from constraints
  select(type, eqn_sys, equation):

  # Convert equations to expressions
  map(x -> lhs(x) - rhs(x), %):

  # Simplify
  expand(%):

  # Remove 0 so we don't divide by 0
  % minus {0}:

for expr in % do

  # Divide each expression in new_eqn_sys by this expression
  map(x -> x/expr, new_eqn_sys):

  # Simplify
  radnormal(%):

  # Check if type is a constant
  map(type, %, constant):

  # Check that ANY are constants
  `or`(op(%)):

  # If not (ie `expr` isn't a constant multiple of any other) add to
  # new_eqn_sys
  if not % then
    new_eqn_sys := new_eqn_sys union {expr}:
  end if:
end do:

  # Determine solving variables
  SolvingVars(eqn_sys, const, ':-ind_vars'=ind_vars):

  # Intersect with ranking
  new_ranking := % intersect indets(ranking):

# Check if this would be better solved in polar coordinates
if numelems(new_eqn_sys) = 2 # Do we have two equations?
  and numelems(new_ranking) = 2 # Do we only wave one new_ranking set?
  and [op(new_ranking[1])] = [op(new_ranking[2])] # Do they have the same functional dependence?
  and convert(LinearAlgebra:-GenerateMatrix(
    Diff2Diffop(new_eqn_sys)[1], new_ranking)[2],list) <> [0,0] then # Are the equations nonlinear in these variables?

  # Split equations and constraints
  constraints := selectremove(type, eqn_sys, equation)[2]:
  pde_equations := map(x -> x = 0, new_eqn_sys):

  # We have two, nonlinear equations; convert to polar form

  # Convert to a set since order will matter
  orig_vars := convert(new_ranking,list):
  params := op(orig_vars[1]):

  # Define new variables
  _A := UniqueName():
  _phi := UniqueName():

  # Replace variables with new variables
  eval(pde_equations,
    {op(0,orig_vars[1]) = unapply(_A(params)*cos(_phi(params)),params),
     op(0,orig_vars[2]) = unapply(_A(params)*sin(_phi(params)),params)}):

  # Add equations together
  %[1] + I*%[2]:
  new_eq := lhs(%)-rhs(%):

  # Determine if denominator could ever be zero
  factor(%):
  denom_eq := denom(%) = 0:
  # Check that the denom equation is possible (and not, for instance
  # 1=0), and run a consistency check
  coulditbe(%) and PDEtools:-ConsistencyTest({%} union constraints):
  # If so, take this as a possible solution
  if % then
    # Solve equation
    {solve(denom_eq, {_A(params),_phi(params)})}:

    # Undo replacement
    sols := map(x -> {orig_vars[1] = eval(_A(params)*cos(_phi(params)),x),
      orig_vars[2] = eval(_A(params)*sin(_phi(params)),x)},%):
  end if:

  # Consider case with denominator non-zero
  new_eq*lhs(denom_eq):
  expand(%):

  # Convert to complex exponentials
  convert(%,exp):
  convert(%,diff):

  # Expand so that exp(-I*phi) becomes exp(I*phi)^{-1}
  expand(%):

unassign('X'):
  # Replace exp with X
  eval(%, {exp(I*_phi(params))=X}):
  # Divide out lowest common power of exp
  %/X^(eval(ldegree(%,X),FAIL=0));
  radsimp(%):
  # Undo replacement
  eval(%, {X=exp(I*_phi(params))}):

  # Separate real and imaginary parts
  solving_system := [evalc(Re(%))=0,evalc(Im(%))=0]:

new_sols := {}:

  # Create list of variables for which to solve
  eval({_A(params),_phi(params)}):

  # Solve system
  sols := sols union {pdsolve(solving_system, %, 'ivars'=ind_vars)}:

  # Replace constant of integration for phi with arctan (so that the
  # constant doesn't appear inside a trig function)
  for sol in % do

temp_sol:=sol:

    # Hack-y: generate a random constant to add to each equation so we can
    # guarantee every equation is a sum and not a monomial; add two
    # constants together in case the equation was 0=0, we need to to still
    # get a sum
    UniqueName():
    unique_const := % + UniqueName():

    eval(_phi(params), sol):

    # Make sure we have a sum of terms
    expand(%):
    % + unique_const:

    select(type, %, Or(function, name)):

    indets(%, assignable):

    % minus existing_names minus {_A(params),_phi(params)}:
    % minus indets(unique_const):

    for repl_name in % do
      if type(repl_name, function) then
        op(repl_name):
        op(0,repl_name) = unapply(arctan(UniqueName()(%)),%):
      else
        repl_name = arctan(UniqueName()):
      end if:
      temp_sol := map(x->lhs(x)=eval(rhs(x),%),sol):
    end do:

    new_sols := {op(new_sols), temp_sol}:

  end do:

  expand(new_sols):

  # Undo replacement
  sols := map(x -> {orig_vars[1] = eval(_A(params)*cos(_phi(params)),x),
    orig_vars[2] = eval(_A(params)*sin(_phi(params)),x)},%):

  # Expand (to get rid of arctan)
  sols := expand(sols):

  # Only keep those constraints that could be satisifed
  sols := map(CheckConstraints, sols, constraints):

  if sols = {} then
    # Try solving like regular. Note: we didn't do this first as Maple
    # will often hang on a nonlinear equation like this. Therefore, we
    # try the faster polar method first
    {pdsolve(pde_equations, new_ranking, ivars=ind_vars)}:

    # Only keep those constraints that could be satisifed
    # Note: we didn't include constraints in pdsolve as they could be
    # non-polynomial in the solving variables (eg sin(f(t)) <> 0)
    sols := map(CheckConstraints, %, constraints):

  end if:

elif ranking = [] then

  unassign('new_ranking'):
  # Create list of variables for which to solve
  SolvingVars(eqn_sys, const, ':-ind_vars'=ind_vars):

  new_ranking := [%]:

  # Solve for all unknown variables
  {pdsolve(eqn_sys, new_ranking, ':-ivars'=ind_vars)}:

else

  # Solve equations

  {pdsolve(eqn_sys, ranking, ':-ivars'=ind_vars)}:

end if:


  # Remove partially solved solutions
  remove(type, %, list):

  # Convert RootOf to radical
  convert(%, radical):
  convert(%, expsincos):

  # Occassionally pdsolve will generate arctan(sin(...)/cos(...)) in the solution; try
  # converting the argument of arctan to tangent and cancelling inverses
  evalindets(%, specfunc(arctan), convert,tan):
  SolveTools:-CancelInverses(%);
  simplify(%):

#  remove(has, %, I):
  remove(has, %, int):
  remove(has, %, RootOf):

  sols := value(%):

  # Determine new variables
  unassign('expr'):
  indets(sols, 'assignable') minus existing_names:

  evalindets(%, function, x -> op(0,x)):

  for expr in % do
    sols := eval(sols, expr = UniqueName()):

  end do:

  return sols:

end proc:

SolveICs := proc(
  ICs::Or(set,list)(equation),
  sols::Or(set,list)(assignable=algebraic),
  const::Or(set,list)(name):={},
  {ind_vars::Or(set,list)({name,function}):={}},
  $)
description "Solve a set/list of intial conditions ICs applied to a set/list of solutions. 'const' is a set/list of variables that will not be solved for. If `ind_vars` is supplied, only functions of these variables will be considered.":

  # Apply initial conditions (don't use exclude_vars)
  EvalSolutions(ICs, sols):

  # If any of the equations cannot be solved, return
  if map(coulditbe, %) <> {true} then
    return {{}}:
  end if:

  # Solve initial conditions
  SolveEquations(%, ':-const'=const, ':-ind_vars'=ind_vars):

  # Only keep solutions that are fully solved
  remove(has, %, int):

  return %:

end proc:

SolvingVars := proc(
  eqn_sys::Or(set,list)(Or(`=`,`<>`)),
  const::Or(set,list)(name):={},
  {ind_vars::Or(set,list)({name,function}):={}},
  $)
local
  indeps:
description "Return a set of variables which appear in set of equations or inequations `eqn_sys` but are not in set `const`. If `ind_vars` is passed, only functions depending on a subset of `ind_var` will be returned.":

  # Convert D to Diff, otherwise indets/assignable won't recognize f(t)
  # in D(f)(t)
  convert(eqn_sys,Diff):

  # Create list of variables for which to solve
  indets(%,`assignable`) minus const minus ind_vars:

  # Remove independent variables (converting to Diff) may have
  # introduced spurious variables. Eg: D(f)(0) converts to
  # eval(diff(f(t1),t1),t1=0), so t1 would be spurious
  % minus IndependentVars(%):

  # Remove duplicate functions
  DedupNamesFuncs(%, ':-ind_vars'=ind_vars):

  # Convert to set: there is no logical ordering to these, so a list
  # wouldn't make sense)
  convert(%, set):

  return %:
end proc:

ToAnalytic := proc(
  expr,
  var,
  $)
local
  s,
  x:
description "convert expression from function of (var) to analytic representation of (var)":

  # Replace var with dummy variable
  subs(var=s, expr):
  # (Inverse) Hilbert transform dummy variable back to var
  inttrans:-invhilbert(%, s, x):
  # Combine trig functions
  combine(%,trig):
  # Undo substitution
  subs(x=var,%):
  # Create analytic representation
  I*% + expr:
  # Convert to complex exponentials
  convert(%,exp):
  # Expand and simplify
  expand(%):
  simplify(%):

  return %:
end proc:

TrigCoeff := proc (
  expr::algebraic,
  basis::function,
  exclude_vars::Or(set,list),
  $)
description "Calculate the coefficient of `basis` trigonometric function in expr after combining trigonometric functions in `expr`. `exclude_vars` will be excluded from trigonometric combination: `exclude_vars` should contain the argument of `basis` function.":

  # Expand first (makes combine/trig quicker)
  expand(expr):

  # Combine powers of sin/cos into multiple-angles form
  combine(%, 'trig'):

  # Expand trig so that any sin(A+B) terms are split
  TrigExpand(%, exclude_vars):

  # Calculate coefficient of `basis`
  coeff(%, basis):

  return %:

end proc:

TrigExpand := proc(
  expr,
  argument_var::Or({name,function},Or(set,list)({name,function})),
  $)
local
  frozen_expr,
  elem:
description "Expand trig functions in `expr` without expanding multiples of `argument_var`":

  # Combine trig so that sin(n*(x+y)) becomes sin(n*x + n*y)
  frozen_expr := evalindets(expr, specfunc({sin,cos}), x -> op(0,x)(expand(op(x))));

  for elem in argument_var do
    # Freeze n*elem
    frozen_expr := subsindets(frozen_expr,
      And(specop(`*`),anyop(integer,identical(elem))), x->freeze(x));
  end do:

  # Expand trig so that any sin(A + B) terms are split
  expand(%):

  thaw(%);

  return %:

end proc:

UnapplyFuncs := proc(
  functions::Or(set,list)({name,function}=algebraic),
  $)
local
  fcn,
  arguments,
  output:
description "Given a set of equalities with names or functions on the LHS and algebraic functionsessions on the RHS, unapply any functions on the LHS":

  output := []:
  for fcn in functions do

    if lhs(fcn) = rhs(fcn) then
      # If LHS = RHS, no point assigning; doing so for functions
      # introduces infinite recursion problems
      next:
    end if:

    # Check if LHS(fcn) is a function (and its arguments are names, not
    # constants; ie, it's not an initial condition)
    if type(lhs(fcn), function) and not type({op(lhs(fcn))},set(constant)) then
      # Get arguments of fcn
      [op(lhs(fcn))]:

      # Define new names for arguments (which could themselves be functions)
      arguments := map(x -> x = UniqueName(), %):

      # Replace arguments with new names
      subs(arguments, rhs(fcn)):

      # Unapply function
      'op(0,lhs(fcn))' = unapply(%, map(rhs, arguments)):

      # Store results
      output := [op(output), %]:
    else
      # Otherwise, just append to ouput
      output := [op(output), fcn]:
    end if:
  end do:

  # If `functions` was a set, convert output to a set as well
  if type(functions, set) then
    output := convert(output, set):
  end if:

  return output:

end proc:

UniqueName := proc()
description "Return a uniquely generated name of the form _D{$num} where $num is increased by on after each call.":

  if not type(counter, integer) then
    # Initialize counter
    counter := 0:
  end if:

  counter := counter+1:
  return cat(_D,counter):
end proc:

CombineSolutions := proc(
  full_solution::Or({name,function}=algebraic,
    {set,list}({name,function}=algebraic)),
  partial_solutions::Or(set,list)({name,function}=algebraic),
  {exclude_vars::Or(set,list):={}},
  $)
description "Use `partial_solutions` to evaluate `full_solution`, possibly excluding the evaluation of `exclude_vars`":

  unassign('j'):

  # Apply `partial_solutions` to `full_solution`
  EvalSolutions(full_solution, partial_solutions, exclude_vars):

  # Calculate excluded solutions
  select((expr,arg) -> member(lhs(expr),arg), partial_solutions,
    exclude_vars):

  return '%%', '%':

end proc:

CompleteSystem := proc(
  eqns::{set,list}(equation),
  expansion_order::integer,
  {scales::name:={},
   scales_order::integer:=expansion_order-leading_order_eqs+1,
   ind_vars::Or(set,list)({name,function}):={},
   expansion_var::name:=epsilon,
   leading_order_eqs::integer:=0,
   const::Or(set,list)(name):={},
   constraints::Or(set,list)({`=`({name,function},constant),`<>`({name,function},constant),
     `<`({name,function},constant),`>`({name,function},constant),
     `<=`({name,function},constant),`>=`({name,function},constant)}):={},
   ICs::Or(set,list)(equation):={},
   exclude_vars::Or(set,list)({function,name}):={},
   assumptions::Or(set,list)(Or(name::property,And(relation,anyop(name,algebraic)))):={},
   extract_type::{"trig","innerProduct"}:="trig",
   solve_eqns::{range,int}:={},
   leading_order_terms::Or(set,list)(list):={},
   secular_allowed::Or(set,list)(function):={},
   restrict_scales::truefalse:=false},
  $)
local
  pdes,
  IndDefs,
  DepDefs,
  sols,
  complementary_scales,
  split_ICs,
  new_ind_vars:
description "Wrapper function for GenEqns, SolveSystem, and VerifyAndCombine. Relevant parameters are documented in those functions.  The only new option is solve_eqs: if eqns is a list and solve_eqns is an integer or range of integers, then only those indexed equations in eqns will be solved (all will be verified though). This is useful if some of the equations can be solved exactly and are used to constrain trail_functions.":

  # Generate order-by-order expansion of equations
  pdes, IndDefs, DepDefs, new_ind_vars, complementary_scales, split_ICs :=
    MultipleScales:-GenEqns(eqns, expansion_order, ':-scales'=scales,
    ':-scales_order'=scales_order, ':-ICs'=ICs,
    ':-leading_order_eqs'=leading_order_eqs, ':-ind_vars'=ind_vars,
    ':-leading_order_terms'=leading_order_terms,
    ':-restrict_scales'=restrict_scales):

  if type(eqns, list) and solve_eqns <> {} then
    # Only pass equations to be solved
    pdes := map(x -> [op(solve_eqns, x)], pdes):
  end if:

  sols := MultipleScales:-SolveSystem(pdes, IndName(IndDefs),
    ':-const'=const, ':-constraints'=constraints,
    ':-exclude_vars'=exclude_vars, ':-ICs'=split_ICs,
    ':-assumptions'=assumptions, ':-extract_type'=extract_type,
    ':-ind_vars'=new_ind_vars, ':-secular_allowed'=secular_allowed,
    ':-complementary_scales'=complementary_scales):

  # Combine solutions and verify equations
  map(x -> MultipleScales:-VerifyAndCombine(eqns, x, expansion_order,
    ':-IndDefs'=IndDefs, ':-DepDefs'=DepDefs,
    ':-exclude_vars'=exclude_vars, ':-assumptions'=assumptions),
    sols):

  return %:

end proc:

GenEqns := proc(
  eqns::{set,list}(equation),
  expansion_order::integer,
  {scales::name:={},
   scales_order::integer:=expansion_order-leading_order_eqs+1,
   ind_vars::Or(set,list)({name,function}):={},
   expansion_var::name:=epsilon,
   ICs::Or(set,list)(equation):={},
   leading_order_eqs::integer:=0,
   leading_order_terms::Or(set,list)([name,nonnegint]):={},
   restrict_scales::truefalse:=false},
  $)

local
  leading_order,
  j,
  pdes,
  new_ind_vars,
  replaced_ind_vars,
  ind_scales,
  ind_leading_order,
  ind_trailing_order,
  IndDef,
  IndFuncExpnd,
  IndNameExpnd,
  ind_ICs,
  DepDef,
  dep_var,
  complementary_scales,
  split_ICs:

description "Given a set/list of differential equations `eqns`, generation a list of sets/lists of equations at each order in a expansion parameter `expansion_var`. If a name `scales` is provided, this variables will be expanded in a perturbation series; otherwise, all independent variables in `eqns` will be expanded in a perturbation series. If `ind_vars` is supplied, only functions of these variables will be considered. If a set/list of equations/inequations `ICs` is provided, they will also be split order-by-order in the expansion parameter. If a set/list of lists `leading_order_terms` of the form {[name,nonnegint],...} is provided, then the expansion of each variable appearing as a name in `leading_order_terms` will have an expansion starting at its corresponding nonnegint value instead of at 0.  Value `name` for functions in leading_order_terms` can be either the function name (eg f) or the function with dependence (eg f(x,t)). If `restrict_scales` is true (default: false) then complementary (homogeneous) solutions for higher-order dependent variables will have depend on fewer indpenedent variable scales. For example, with expansion_order=2, a dependent variable y(t) would be expanded as y__0(t__0,t__1,t__2) + y__1(t__0,t__1) + y__2(t__0). This prevents higher-order terms from having undetermined homogeneous solutions.":

  # Check if scales was passed
  if scales = {} then
    # Get a list of dependent variables
    ind_scales := IndependentVars(eqns):
    if numelems(%) = 0 then
      error "Could not determine independent variable. Must pass `scales` denoting independent variable.":
    elif numelems(%) > 1 then
      error "More than one independent variable. Must pass `scales` denoting independent variable to expand.":
    end if:
  else
    # Otherwise, just use scales
    ind_scales := scales:
  end if:

  # Determine starting order for independent variable
  ind_leading_order := 0:
  # Check if leading_order_terms has an entry for scales
  select(elem -> op(1,elem)=scales, leading_order_terms):

  if numelems(%) > 1 then
    error cat("More than one entry for ",scales," in leading_order_terms dictionary."):
  elif numelems(%) = 1 then
    ind_leading_order := op([1,2],%):
  end if:

  # Define expansion for independent variable
  IndDef := IndGen(ind_scales, scales_order+ind_leading_order,
    ':-expansion_var'=expansion_var,
    ':-leading_order_terms'=ind_leading_order):
  ind_ICs := IndICs(ind_scales, scales_order+ind_leading_order,
    ':-expansion_var'=expansion_var,
    ':-leading_order_terms'=ind_leading_order):

  # Define replacement for independent variable with its expansion
  unassign('j'):
  ind_trailing_order := ind_leading_order+scales_order-1:
  IndFuncExpnd := map(x -> x = seq(cat(x,__,j)(x),
    j=ind_leading_order..ind_trailing_order ), ind_scales):
  IndNameExpnd := map(x -> x = seq(cat(x,__,j),
    j=ind_leading_order..ind_trailing_order), ind_scales):

  # Generate list of all independent variables
  new_ind_vars := {op(scales), op(ind_vars)}:

  # Generate list of dependent variables
  DepDef := {}:
  complementary_scales := {}:

  DependentVars(eqns, ':-ind_vars'=%):

  for dep_var in % do
    # Determine starting order
    leading_order := 0:
    # Check if leading_order_terms has an entry for scales
    select(elem ->op(1,elem)=dep_var or op(1,elem)=op(0,dep_var), leading_order_terms):
    if numelems(%) > 1 then
      error cat("More than one entry for ",dep_var," in leading_order_terms dictionary."):
    elif numelems(%) = 1 then
      leading_order := op([1,2],%):
    end if:

    # Replace independent scales with their # expansions
    replaced_ind_vars := op(eval(dep_var,IndFuncExpnd)):
    # Add dep_var to set
    DepDef := DepDef union {dep_var =
      add(cat(op(0,dep_var),__,j)(replaced_ind_vars)
      *expansion_var^j, j=leading_order..expansion_order)}:

    if restrict_scales then
      # Restrict operands of dep_var
      {seq([cat(op(0,dep_var),__,j)(replaced_ind_vars), [op(eval(dep_var,
	lhs(IndFuncExpnd) =
	(op([ind_leading_order+1..min(ind_trailing_order+1,
        expansion_order-j+1)], [rhs(IndFuncExpnd)]))))]
        ], j=leading_order..expansion_order)}:

      complementary_scales := complementary_scales union %:
    end if:

  end do:

  # Convert independent functions to names; this makes solving the systems easier
  complementary_scales := eval['recurse'](complementary_scales,
    IndFuncToName(IndDef)):

  # Make substitutions (leave independent variables unevaluated)
  EvalSolutions(eqns, [op(DepDef), op(IndDef)], IndFunc(IndDef)):

  # Convert independent functions to names; this makes solving the systems easier
  eval['recurse'](%, IndFuncToName(IndDef)):

  # Remove terms above order expansion_order
  map2(map,TruncateExpr, % ,expansion_order,
    ':-expansion_var'=expansion_var):

  # Collect terms order-by-order
  pdes := map[5](map2, map, coeff, %, expansion_var, [seq(j,j=leading_order_eqs..expansion_order)]):

  # Replace scales with it their expansions in list of all independent variables
  new_ind_vars := eval(new_ind_vars, IndNameExpnd):

  # Expand ICs
  split_ICs := []:
  if ICs <> {} then
     # Make substitutions (leave independent variables unevaluated)
    EvalSolutions(ICs, [op(DepDef), op(IndDef)], IndFunc(IndDef)):

    # Convert independent functions to names; this makes solving the
    # systems easier
    eval['recurse'](%, ind_ICs):
    value(%):

    # Remove terms above order expansion_order
    map2(map,TruncateExpr, % ,expansion_order,
      ':-expansion_var'=expansion_var):

    # Collect terms order-by-order
    map[5](map2, map, coeff, %, expansion_var,
      [seq(j,j=leading_order_eqs..expansion_order)]):

    # Order of equations doesn't matter, so combine them all into one
    # set
    map(convert, %, set):
    split_ICs := `union`(op(%)):
  end if:

  return pdes, IndDef, DepDef, new_ind_vars, complementary_scales,
    split_ICs:

end proc:

IndICs := proc(
  scales::(name),
  scales_order::integer,
  {expansion_var::name:=epsilon,
   leading_order_terms::nonnegint:=0},
  $)
local
  ind_scales:

description "Given a name `scales` and an integer `scales_order`, return a lists of `scales_order`-number of functions for extraction initial conditions by extracting the appropriate perturbation order. For example, if scales=t and scales_order=3, then return will be [t__0(t) = t->%coeff(t,epsilon,0), t__1(t) = t->%coeff(t,epsilon,1), t__2(t) = t->%ceff(t,epsilon,2). This way, the independent variable can be broken into its constituent expansions: t__1(A*epsilon + B) = A, t__0(epsilon^2)=0, etc. This is useful for applying initial conditions, where---for instance---t=0 should imply 0=t__0=t__1=t__2=... . If a non-negative integer `leading_order_terms` is given, the expansion starts at that integer; otherwise, it starts at 0.":

  # Check if scales was passed
  if scales = {} then
    # Get a list of dependent variables
    ind_scales := IndependentVars(eqns):
  else
    # Otherwise, just use scales
    ind_scales := scales:
  end if:

  unassign('j'):
  # Define multiple scales
  map(x -> [seq(cat(x,__,j) = unapply(%coeff(x,expansion_var,j),x),
    j=leading_order_terms..scales_order-1)], ind_scales):

  return %:

end proc:

IndFunc := proc(
  IndDef::Or(set,list)(function=algebraic),
  $)
description "Given a set/list of equations (function=algebraic), returns a set/list of the functions.":

  map(lhs, IndDef):

  return %:

end proc:

IndFuncToName := proc(
  IndDef::Or(set,list)(function=algebraic),
  $)

description "Given a set/list of equations (function=algebraic), returns a set/list of equations equating the functions to the names of the functions.":

  map(x -> op(IndFunc({x})) = op(IndName({x})), IndDef):

  return %:

end proc:

IndName := proc(
  IndDef::Or(set,list)(function=algebraic),
  $)

description "Given a set/list of equations (function=algebraic), returns a set/list of the names of the functions.":

  IndFunc(IndDef):

  map2(op, 0, %):

  return %:

end proc:

IndGen := proc(
  scales::(name),
  scales_order::integer,
  {expansion_var::name:=epsilon,
   leading_order_terms::nonnegint:=0},
  $)
local
  ind_scales:

description "Given a name `scales` and an integer `scales_order`, return a lists of `scales_order`-number of definitions for perturbation scales in `expansion_var`. For example, if scales=t and scales_order=3, then return will be [t__0(t) = t, t__1(t) = epsilon*t, t__2(t) = epsilon^2*t]. If a non-negative integer `leading_order_terms` is given, the expansion starts at that integer; otherwise, it starts at 0.":

  # Check if scales was passed
  if scales = {} then
    # Get a list of dependent variables
    ind_scales := IndependentVars(eqns):
  else
    # Otherwise, just use scales
    ind_scales := scales:
  end if:

  unassign('j'):
  # Define multiple scales
  map(x -> [seq(cat(x,__,j)(x) = x*expansion_var^j,
    j=leading_order_terms..scales_order-1)], ind_scales):

  return %:

end proc:

IndNameToFunc := proc(
  IndDef::Or(set,list)(function=algebraic),
  $)

description "Given a set/list of equations (function=algebraic), returns a set/list of equations equating the names of the functions to the functions.":

  IndFuncToName(IndDef):

  map(x -> rhs(x) = lhs(x), %):

  return %:


end proc:

SolveSystem := proc(
  pdes::list({set,list}(Or(`=`, `<>`))),
  mult_scales::list(name),
  {const::Or(set,list)(name):={},
   ind_vars::Or(set,list)({name,function}):={},
   exclude_vars::Or(set,list)({function,name}):={},
   constraints::Or(set,list)({`=`({name,function},constant),`<>`({name,function},constant),
     `<`({name,function},constant),`>`({name,function},constant),
     `<=`({name,function},constant),`>=`({name,function},constant)}):={},
   ICs::Or(set,list)(equation):={},
   assumptions::Or(set,list)(Or(name::property,And(relation,anyop(name,algebraic)))):={},
   extract_type::{"trig","innerProduct"}:="trig",
   secular_allowed::Or(set,list)(function):={},
   complementary_scales::{set,list}([{name,function},list(name)]):={}},
  $)
local
  sols,
  sol,
  new_sols,
  elem,
  n:
description "Solve a system of partial differential equations (PDEs) perturbatively. `pdes` must be a list, with the `n`th element being a set/list of `n`-th order PDEs. `solving_order` denotes which order in the perturbative parameter `PDEs` satisfies. `mult_scales` is a list of names of independent variables; secular terms in these variables will be avoided (in order) unless the term becoming secular appears in `allowed_secular`. `const` is a set of constants that should not be solved for (arbitrary constants of the problem).  `exclude_vars` are functions/variables whose values should not be substituted into the answer (though derivative of these quantities will still be evaluated).  Additionally, `constraints` are any additional equation/inequalities of the form 'function/symbol =/</>/<=/>=/<> constant'. `ICs` are any initial conditions to be imposed; ICs should already be split by perturbation order (ie the perturbation parameter should not appear in ICs). `assumptions` are any assumed properties to be applied when choosing which of multiple solutions to keep (eg, some parameters are real). `extract_type` is passed to ExtractCoeffs; see ExtractCoeffs for a description. If `ind_vars` is supplied, only functions of these variables will be considered. If a set/list of lists `complementary_scales` of the form [{name,function},[names]] is supplied, then any function appearing in the first element (either a name or function) will only have homogenous terms depending on the variables `names`. This is useful for restricting higher-order terms from having too many free parameters.":

  # Create solution set
  sols := {{}}:

  # Solve equations order-by-order
  unassign('n'):
  for n from 1 to numelems(pdes) do

    # Create empty set to hold new solutions
    new_sols := {}:

    for sol in sols do
      # Solve this order
      SolveOrder(pdes[n], n, mult_scales, ':-const'=const,
        ':-prev_sols'=sol, ':-exclude_vars'=exclude_vars,
        ':-constraints'=constraints, ':-extract_type'=extract_type,
        ':-ind_vars'=ind_vars, ':-secular_allowed'=secular_allowed,
        ':-complementary_scales'=complementary_scales):

      # If no solution, go to next element
      if % = {{}} then next end if:

      # Evaluate previous solutions with new solutions and append new
      # solutions
      map(sol_set -> EvalSolutions(sol, sol_set, exclude_vars) union
        sol_set, %):

      # Add solutions to current set
      new_sols := new_sols union %:

    end do:

    # Make sure solution set is not empty
    if numelems(new_sols) = 0 then
      error cat(`No solution found at order `, n):
    end if:

    # Replace solutions with newly created solution set
    sols := new_sols:

  end do:

  if numelems(ICs) <> 0 then

    # Create empty set to hold new solutions
    new_sols := {}:

    # Determine all variables to solve for
    indets(sols, assignable(name)) minus ind_vars minus const:

    # Add an empty set first
    [{}, op(%)]:

    # Sometimes, a constant of integration needs to be taken to infinity
    # to solve an initial condition; pdsolve can't figure this out, so
    # manually do it
    for elem in % do

      # One at a time, try taking each element to infinity
      if elem <> {} then
        map2(map,limit, sols, elem=infinity):
      else
        sols:
      end if:

      # If any terms are undefined, try the next term
      if has(%, undefined) or has(%, infinity) then
        next:
      end if:

      for sol in % do
        # Apply initial conditions
        SolveICs(ICs, sol, const, ':-ind_vars'=ind_vars):

        # If no solution, go to next element
        if % = {{}} then next end if:

        # Evaluate previous solutions with new solutions and append new
        # solutions
        map(sol_set -> EvalSolutions(sol, sol_set, exclude_vars) union
          sol_set, %):

        # Add solutions to current set
        new_sols := new_sols union %:

      end do:

      # If we found a solution, continue
      if numelems(new_sols) <> 0 then
        break:
      end if:

    end do:

    # Make sure solution set is not empty
    if numelems(new_sols) = 0 then
      error "No solution found to initial conditions":
    end if:

    # Replace solutions with newly created solution set
    sols := new_sols:

  end if:

  # `sols` now contains many auto-generated terms. Only return the
  # solutions to terms originally appearing in PDEs
  indets(pdes, assignable) minus ind_vars minus const:

  # Evaluate terms using solutions
  map[3](map, (x,y) -> x = eval['recurse'](x, y), %, sols):

  # Return solutions
  return %:

end proc:

TruncateExpr := proc(
  expr::algebraic,
  trunc_order::integer,
  {expansion_var::name:=epsilon},
  $)
description "Truncate `expr` to order `trunc_order` in variable `expansion_var`":

  # Remove terms above order `trunc_order`
  convert(taylor(expr,expansion_var=0,trunc_order+1),polynom):

  return %:
end proc:

UnitTests := proc()

local
  pdes,
  const,
  ind_vars,
  constraints,
  exclude_vars,
  mult_scales,
  ICs,
  PDE,
  solving_order,
  assumptions,
  complementary_scales,
  extract_type,
  sols:

description "Run unit tests.":

# Test TruncateExpr
CodeTools:-Test(TruncateExpr(epsilon^3*F(x) + epsilon*G(x),
  2), epsilon*G(x) , label="TruncateExpr Test 1"):
CodeTools:-Test(TruncateExpr(a + g*x + g^2 + F(g^5) - g^4*y, 3,
  'expansion_var'=g), a + g*x + g^2 + F(0), label="TruncateExpr Test 2"):

# Test EvalSolutions
CodeTools:-Test(EvalSolutions({f(x)*g(x)},{f(x) = x^2}), {x^2*g(x)},
  label="EvalSolutions Test 1"):
CodeTools:-Test(EvalSolutions({f(g(x))+diff(g(x),x)},{g(x) = A*x^3},
  {g(x)}), {f(g(x)) + 3*A*x^2}, label="EvalSolutions Test 2"):
CodeTools:-Test(EvalSolutions({Eval(cos(f(x))*D(f)(x),f(x)=0)}, {}),
  {D(f)(x)}, label="EvalSolutions Test 3"):
CodeTools:-Test(EvalSolutions({D[1,1](f)(t) + f(t)},{f(t) = t^2},
  {f(t)}), {2 + f(t)}, label="EvalSolutions Test 4"):

# Test TrigExpand
CodeTools:-Test(TrigExpand(cos(2*x + y), {x}),
  cos(2*x)*cos(y) - sin(2*x)*sin(y),
  label="TrigExpand Test 1"):
CodeTools:-Test(TrigExpand(cos(x-y), {y}), cos(y)*cos(x)+sin(y)*sin(x),
  label="TrigExpand Test 2"):
CodeTools:-Test(TrigExpand(cos(x-2*y), {y}),
  cos(x)*cos(2*y)+sin(x)*sin(2*y), label="TrigExpand Test 3"):
CodeTools:-Test(TrigExpand(p*cos(-tau(t)+psi),
  {tau(t)}), p*cos(tau(t))*cos(psi)+p*sin(tau(t))*sin(psi),
  label="TrigExpand Test 3"):

# Test BasisCoeff
CodeTools:-Test(BasisCoeff(A*sin(x), sin(x)), A,
  label="BasisCoeff Test 1"):
CodeTools:-Test(BasisCoeff(B*sin(x) + C*cos(y), cos(y)), C,
  label="BasisCoeff Test 2"):
CodeTools:-Test(BasisCoeff(F*sin(x), sin(x), 'limits'=Pi), F,
  label="BasisCoeff Test 3"):
CodeTools:-Test(BasisCoeff(G*cos(z), cos(z), z, 'limits'=Pi), G,
  label="BasisCoeff Test 4"):
CodeTools:-Test(BasisCoeff(H*cos(z)+7, 1, z, 'limits'=Pi), 14,
  label="BasisCoeff Test 5"):

# Test bug fixes
CodeTools:-Test(
  expand(inttrans:-hilbert(sin(a)*sin(t+b), t, s)),
  -sin(s)*sin(a)*sin(b)+cos(s)*sin(a)*cos(b),
  label="inttrans:-hilbert Bug Fix Test"):
CodeTools:-Test(
  PDEtools:-Solve({diff(A(t),t) = A(t)*x}, {A(t)}, ivars={x,t}),
  {A(t)=0},
  label="PDEtools:-Solve Pass-Through Parameters Bug Fix Test"):
CodeTools:-Test(
  PDEtools:-Solve({D(B)(x) = 0, A()* D(B)(x) = 0}, [{A(), B(x)}]),
  {A() = A(), B(x) = _C1},
  label="PDEtools:-Solve Empty Function Bug Fix Test"):
CodeTools:-Test(
  convert(D[2](eta)(t,x)+D[2](phi)(t,x,0), Diff),
  Diff(eta(t,x),x) + Diff(phi(t,x,0),x),
  label="Convert/Diff and Eval Bug Fix Test"):

# Test SolveSystem using 1st order Duffing oscillator
pdes := [{y__0(t__0,t__1)+Diff(Diff(y__0(t__0,t__1),t__0),t__0) = 0},
{y__0(t__0,t__1)^3*nu+y__1(t__0,t__1)+2*Diff(Diff(y__0(t__0,t__1),t__0),t__1)+Diff(Diff
(y__1(t__0,t__1),t__0),t__0) = 0}]:
mult_scales := [t__0,t__1]:

const := {nu}:
constraints := {y__0(t__0,t__1) <> 0}:
ind_vars := {t__0,t__1}:
ICs := {D[2](y__0)(0,0)+D[1](y__1)(0,0) = 0, y__0(0,0) = 1, y__1(0,0) = 0,
  D[1](y__0)(0,0) = 0}:

assumptions := {nu::positive}:
complementary_scales := {[y__0(t__0,t__1), [t__0, t__1]],
  [y__1(t__0,t__1), [t__0]]}:

CodeTools:-Test(expand(SolveSystem(pdes, mult_scales,
  'const'=const, 'constraints'=constraints, 'ind_vars'=ind_vars,
  'ICs'=ICs, 'complementary_scales'=complementary_scales,
  'assumptions'=assumptions)),
{{y__0(t__0,t__1)=-sin(3/8*nu*t__1)*sin(t__0)+cos(3/8*nu*t__1)*cos(t__0),
y__1(t__0,t__1)=-1/32*cos(t__0)*nu+1/8*nu*cos(9/8*nu*t__1)*cos(t__0)^3
-3/32*nu*cos(t__0)*cos(9/8*nu*t__1)-1/8*nu*sin(9/8*nu*t__1)*sin(t__0)*cos(t__0)^2
+1/32*nu*sin(t__0)*sin(9/8*nu*t__1)}}, label="SolveSystem 1st-Order Duffing Test"):

# Test VerifySolutions using 3rd order Duffing equation
PDE := {diff(diff(Y(t),t),t)+Y(t)+epsilon*nu*Y(t)^3 = 0}:
solving_order := 3:
sols := [Y(t) =
a*cos(tau(t__0(t),t__1(t),t__2(t),t__3(t)))+(-1/32*a^3*nu*cos(tau(t__0(t),t__1(t),
t__2(t),t__3(t)))+1/32*a^3*nu*cos(3*tau(t__0(t),t__1(t),t__2(t),t__3(t))))*epsilon+(23/1024*a^5*nu^2*cos(tau(t__0(t),t__1(t),t__2(t),t__3(t)))-3/
128*a^5*nu^2*cos(3*tau(t__0(t),t__1(t),t__2(t),t__3(t)))+1/1024*nu^2*a^5*cos(5*tau(t__0(t),t__1(t),t__2(t),t__3(t))))*epsilon^2+(-547/32768*a^7*
nu^3*cos(tau(t__0(t),t__1(t),t__2(t),t__3(t)))+297/16384*a^7*nu^3*cos(3*tau(t__0(t),t__1(t),t__2(t),t__3(t)))-3/2048*nu^3*a^7*cos(5*tau(t__0(t)
,t__1(t),t__2(t),t__3(t)))+1/32768*nu^3*a^7*cos(7*tau(t__0(t),t__1(t),t__2(t),t__3(t))))*epsilon^3,
tau(t__0(t),t__1(t),t__2(t),t__3(t)) = t__0(t
)+3/8*nu*a^2*t__1(t)-21/256*a^4*nu^2*t__2(t)+81/2048*a^6*nu^3*t__3(t)+_C1,
t__0(t) = t, t__1(t) = epsilon*t, t__2(t) = epsilon^2*t, t__3(t) =
epsilon
^3*t]:

CodeTools:-Test(VerifySolution(PDE, solving_order,
sols), NULL, testnoerror, label="VerifySolution Test"):

# Test SolveSystem with order-1 unforced Stokes Wave
pdes := [[(-phi__k1S__1(t__0)*Pi+D(eta__k1S__1)(t__0)*Pi)/Pi = 0,
  (-phi__k1C__1(t__0)*Pi+D(eta__k1C__1)(t__0)*Pi)/Pi = 0, 0 = 0,
  (D(phi__k1S__1)(t__0)
  *Pi+eta__k1S__1(t__0)*Pi)/Pi = 0,
  (D(phi__k1C__1)(t__0)*Pi+eta__k1C__1(t__0)*Pi)/Pi = 0,
  D(phi__k0__1)(t__0) = 0, phi__k0__1(t__0) = 0]]:

mult_scales := [t__0]:

ind_vars := {t__0, x, z}:

constraints := {eta__k1S__1(t__0) <> 0}:
ICs := {eta__k1C__1(0) = 0, eta__k1S__1(0) = 0,
  D(eta__k1C__1)(0) = 1, D(eta__k1S__1)(0) = 1}:

complementary_scales := {[eta__k1C__1(t__0), [t__0]],
  [eta__k1S__1(t__0), [t__0]], [phi__k0__1( t__0), [t__0]],
  [phi__k1C__1(t__0), [t__0]], [phi__k1S__1(t__0), [t__0]]}:

CodeTools:-Test(expand(SolveSystem(pdes, mult_scales,
  'constraints' = constraints, 'ICs' = ICs, 'ind_vars' = ind_vars,
  'complementary_scales' = complementary_scales)),
  {{eta__k1C__1(t__0) = sin(t__0), eta__k1S__1(t__0) = sin(t__0),
  phi__k0__1(t__0) = 0, phi__k1C__1(t__0) = cos(t__0), phi__k1S__1(t__0)
  = cos(t__0)}}, label="SolveSystem 1st-Order StokesWaves Test"):

# Test SolveSystem with order-2 unforced Stokes Wave
pdes := [[D[1](eta__k0__1)(t__0,t__1) = 0,
  -phi__k1S__1(t__0,t__1)+D[1](eta__k1S__1)(t__0,t__1) = 0,
  -phi__k1C__1(t__0,t__1)+D[1](eta__k1C__1)(t__0,t__1) = 0 ,
  eta__k0__1(t__0,t__1) = 0,
  2*D[1](phi__k0__1)(t__0,t__1)+eta__k1S__1(t__0,t__1)+D[1](phi__k1S__1)(t__0,t__1)
  = 0, eta__k1C__1(t__0, t__1)+D[1](phi__k1C__1)(t__0,t__1) = 0,
  -D[1](phi__k0__1)(t__0,t__1) = 0, D[1](phi__k0__1)(t__0,t__1) =
  0, phi__k0__1(t__0,t__1) = 0], [D[1](
  eta__k0__2)(t__0,t__1)+D[2](eta__k0__1)(t__0,t__1) = 0,
  -phi__k1S__2(t__0,t__1)+D[2](eta__k1S__1)(t__0,t__1)+D[1](eta__k1S__2)(t__0,t__1)-phi__k0__1(
  t__0,t__1)*eta__k1C__1(t__0,t__1)-eta__k0__1(t__0,t__1)*phi__k1S__1(t__0,t__1)
  = 0, -phi__k1C__2(t__0,t__1)+D[2](eta__k1C__1)(t__0,t__1)+D[1](
  eta__k1C__2)(t__0,t__1)-eta__k0__1(t__0,t__1)*phi__k1C__1(t__0,t__1)+phi__k0__1(t__0,t__1)*eta__k1S__1(t__0,t__1)
  = 0, D[1](eta__k2S__2)(t__0,t__1)-
  phi__k1C__1(t__0,t__1)*eta__k1S__1(t__0,t__1)-phi__k1S__1(t__0,t__1)*eta__k1C__1(t__0,t__1)-2*phi__k2S__2(t__0,t__1)
  = 0,
  D[1](eta__k2C__2)(t__0,t__1)+phi__k1S__1(t__0,t__1)*eta__k1S__1(t__0,t__1)-phi__k1C__1(t__0,t__1)*eta__k1C__1(t__0,t__1)-2*phi__k2C__2(t__0,t__1)
  = 0,
  1/2*phi__k1S__1(t__0,t__1)^2+1/2*phi__k1C__1(t__0,t__1)^2+1/2*phi__k0__1(t__0,t__1)^2+1/2*D[1](phi__k1C__1)(t__0,t__1)*eta__k1C__1(t__0,t__1)+eta__k0__2(t__0,t__1)+1/2*D[1](
  phi__k1S__1)(t__0,t__1)*eta__k1S__1(t__0,t__1) = 0,
  2*D[2](phi__k0__1)(t__0,t__1)+2*D[1](phi__k0__2)(t__0,t__1)+D[1](phi__k1S__1)(t__0,t__1)*
  eta__k0__1(t__0,t__1)-phi__k1C__1(t__0,t__1)*phi__k0__1(t__0,t__1)+eta__k1S__2(t__0,t__1)+D[2](phi__k1S__1)(t__0,t__1)+D[1](phi__k1S__2)(t__0,t__1)
  = 0,
  D[1](phi__k1C__2)(t__0,t__1)+D[1](phi__k1C__1)(t__0,t__1)*eta__k0__1(t__0,t__1)+phi__k1S__1(t__0,t__1)*phi__k0__1(t__0,t__1)+eta__k1C__2(t__0,t__1)+D[2](phi__k1C__1)(t__0,t__1)
  = 0,
  -D[2](phi__k0__1)(t__0,t__1)-D[1](phi__k0__2)(t__0,t__1)+1/2*D[1](phi__k1C__1)(t__0,t__1)*eta__k1S__1(t__0,t__1)+
  1/2*D[1](phi__k1S__1)(t__0,t__1)*eta__k1C__1(t__0,t__1)+D[1](phi__k2S__2)(t__0,t__1)+eta__k2S__2(t__0,t__1)
  = 0, 1/2*D[1](phi__k1C__1)(t__0,t__1)*
  eta__k1C__1(t__0,t__1)+D[1](phi__k2C__2)(t__0,t__1)+eta__k2C__2(t__0,t__1)-1/2*D[1](phi__k1S__1)(t__0,t__1)*eta__k1S__1(t__0,t__1)
  = 0, D[2]( phi__k0__1)(t__0,t__1)+D[1](phi__k0__2)(t__0,t__1) = 0,
  phi__k0__2(t__0,t__1) = 0]]:
mult_scales := [t__0, t__1]:

ind_vars := {t__0, t__1, x, z}:
constraints := {eta__k1S__1 (t__0,t__1) <> 0}:

ICs := {0 = 0, D[2](eta__k1C__1)(0,0)+D[1](eta__k1C__2)(0,0) = 0,
  D[2](eta__k1S__1)(0,0)+D[1]( eta__k1S__2)(0,0) = 0, eta__k1C__1(0,0) =
  0, eta__k1C__2(0,0) = 0, eta__k1S__1(0,0) = 1, eta__k1S__2(0,0) = 0,
  eta__k2C__2(0,0) = 0, eta__k2S__2(0,0) = 0, D[1](eta__k1C__1)(0,0) = 1,
  D[1](eta__k1S__1)(0,0) = 0, D[1](eta__k2C__2)(0,0) = 0,
  D[1](eta__k2S__2)(0,0) = 0}:

complementary_scales := {[eta__k0__1(t__0,t__1), [t__0, t__1]],
  [eta__k0__2(t__0,t__1), [t__0]], [eta__k1C__1(t__0,t__1), [t__0, t__1]],
  [eta__k1C__2(t__0,t__1), [t__0]], [eta__k1S__1(t__0,t__1), [t__0, t__1]],
  [eta__k1S__2(t__0,t__1), [t__0]], [eta__k2C__2(t__0,t__1),
  [t__0]], [eta__k2S__2(t__0,t__1), [t__0]], [phi__k0__1(t__0,t__1),
  [t__0, t__1]], [ phi__k0__2(t__0,t__1), [t__0]],
  [phi__k1C__1(t__0,t__1), [t__0, t__1]], [phi__k1C__2(t__0,t__1),
  [t__0]], [phi__k1S__1(t__0,t__1), [t__0, t__1]],
  [phi__k1S__2(t__0,t__1), [t__0]], [phi__k2C__2(t__0,t__1), [t__0]],
  [phi__k2S__2(t__0,t__1), [t__0]]}:

CodeTools:-Test(expand(SolveSystem(pdes, mult_scales, 'constraints' =
  constraints, 'ICs' = ICs, 'ind_vars' = ind_vars,
  'complementary_scales' = complementary_scales)),
  {{eta__k0__1(t__0,t__1) = 0, eta__k0__2(t__0,t__1) = 0,
  eta__k1C__1(t__0,t__1) = sin(t__0), eta__k1C__2(t__0,t__1) = 0,
  eta__k1S__1(t__0,t__1) = cos(t__0), eta__k1S__2(t__0,t__1) = 0,
  eta__k2C__2(t__0,t__1) = 1/2*cos(2^(1/2)*t__0)-cos(t__0)^2+1/2,
  eta__k2S__2(t__0,t__1) =
  -1/2*2^(1/2)*sin(2^(1/2)*t__0)+sin(t__0)*cos(t__0),
  phi__k0__1(t__0,t__1) = 0, phi__k0__2(t__0,t__1) = 0,
  phi__k1C__1(t__0,t__1) = cos(t__0), phi__k1C__2(t__0,t__1) = 0,
  phi__k1S__1(t__0,t__1) = -sin(t__0), phi__k1S__2(t__0,t__1) = 0,
  phi__k2C__2(t__0,t__1) = -1/4*2^(1/2)*sin(2^(1/2)*t__0),
  phi__k2S__2(t__0,t__1) = -1/2*cos(2^(1/2)*t__0)}},
  label="SolveSystem 2nd-Order StokesWaves Test"):

end proc:

VerifyAndCombine := proc(
  PDEs::Or(set,list)(equation),
  solutions::Or(set,list)(assignable=algebraic),
  e_order::integer,
  {IndDefs::Or(set,list)(function=algebraic):={},
   DepDefs::Or(set,list)(function=algebraic):={},
   exclude_vars::Or(set,list)({function,name}):={},
   assumptions::Or(set,list)(Or(name::property,And(relation,anyop(name,algebraic)))):={}},
  $)
local
  sols,
  gen_sol,
  exclude_sol:
description "Given a set/list of PDEs and a set/list of solutions, verifies that the PDEs are satisifed up to order `e_order`. If PDEs is expressed in terms of independent variables but solutions are expressed in terms of perturbation expansions of independent variales, the definition of the full ind. variables in terms of the expansions should be passed to `IndDefs`; likewise for `DepDefs` for dependent functions. If some variables/functions should be left unevaluated, they should be passed to `exclude_vars`. Finally, any assumptions on variables (eg. positive or real) should be passed to assumptions.":

  # Simplifiy solutions
  map(simplify, solutions) assuming op(assumptions):
  sols := map(combine, %, 'trig'):

  # Combine solutions
  gen_sol, exclude_sol := CombineSolutions(DepDefs, sols,
    ':-exclude_vars'=exclude_vars):

  # Re-add time dependence to exclude_sol (gen_sol should already have it)
  exclude_sol := eval(exclude_sol, IndNameToFunc(IndDefs)):

  # Package solutions together
  sols := [op(gen_sol), op(exclude_sol)]:

  # Re-insert time dependence to ts
  sols := EvalSolutions(sols, IndFuncToName(IndDefs), IndFunc(IndDefs)):

  # Verify solution
  VerifySolution(PDEs, e_order, [op(sols), op(IndDefs)]):

  # Undefine t-dependence (so t[j] are just variables, not functions of
  # t); this makes the output look nicer
  eval['recurse'](sols, IndFuncToName(IndDefs)):

  return %:

end proc:

VerifySolution := proc(
  pdes::{set,list}(Or(`=`, `<>`)),
  solving_order::integer,
  solutions::list(assignable=algebraic),
  {expansion_var::name:=epsilon,
   assumptions::Or(set,list)(Or(`<`,`<=`,`<>`,property,name::property,And(relation,anyop(name,algebraic)))):={}},
  $)
local
  simplified_pdes:
description "Check if the provided `solutions` satisfy the pdes up to order `solving_order` in the perturbative variable `expansion_var` (default `e`). Any `assumptions` are applied while evaluating.":

  # Evaluate pdes using solutions
  simplified_pdes := EvalSolutions(pdes, solutions):

  # Remove terms above order solving_order
  map2(map,TruncateExpr,%, solving_order, ':-expansion_var'=expansion_var):

  # Simplify expressions
  convert(%, exp):
  expand(%):
  simplify(%):

  # Check if equations are TRUE
  map(evalb, %):

  # Convert to set
  convert(%,set):

  if(% <> {true}) then
    error cat(`False solution at `,solving_order,`th order`):
    quit:
  end if:

  return:

end proc:

end module: