From 8c66155128b9836565154ca1e8b98cce72f14cad Mon Sep 17 00:00:00 2001
From: Nathan Fulton <nathan@ibm.com>
Date: Tue, 22 Jun 2021 12:08:44 -0400
Subject: [PATCH 01/14] taylorize interface

---
 .../cs/ls/keymaerax/launcher/KeYmaeraX.scala  | 12 +++-
 .../ls/keymaerax/launcher/TaylorizeMain.scala | 63 +++++++++++++++++++
 2 files changed, 74 insertions(+), 1 deletion(-)
 create mode 100644 keymaerax-webui/src/main/scala/edu/cmu/cs/ls/keymaerax/launcher/TaylorizeMain.scala

diff --git a/keymaerax-webui/src/main/scala/edu/cmu/cs/ls/keymaerax/launcher/KeYmaeraX.scala b/keymaerax-webui/src/main/scala/edu/cmu/cs/ls/keymaerax/launcher/KeYmaeraX.scala
index d29f24c692..d8f7b60a3a 100644
--- a/keymaerax-webui/src/main/scala/edu/cmu/cs/ls/keymaerax/launcher/KeYmaeraX.scala
+++ b/keymaerax-webui/src/main/scala/edu/cmu/cs/ls/keymaerax/launcher/KeYmaeraX.scala
@@ -67,7 +67,8 @@ object KeYmaeraX {
     val CONVERT: String = edu.cmu.cs.ls.keymaerax.cli.KeYmaeraX.Modes.CONVERT
     val SETUP: String = edu.cmu.cs.ls.keymaerax.cli.KeYmaeraX.Modes.SETUP
     val UI: String = "ui"
-    val modes: Set[String] = Set(CODEGEN, CONVERT, MODELPLEX, PROVE, REPL, UI, SETUP)
+    val TAYLORIZE: String = "taylorize"
+    val modes: Set[String] = Set(CODEGEN, CONVERT, MODELPLEX, PROVE, REPL, UI, SETUP, TAYLORIZE)
   }
 
   /** Usage -help information. */
@@ -95,6 +96,13 @@ object KeYmaeraX {
     try {
       //@todo allow multiple passes by filter architecture: -prove bla.key -tactic bla.scal -modelplex -codegen
       options.get('mode) match {
+        case Some(Modes.TAYLORIZE) => {
+          val filename = options.get('file)
+          filename match {
+            case Some(s) => println(TaylorizeMain(s.asInstanceOf[String]))
+            case None => println("FAILED.")
+          }
+        }
         case Some(Modes.CODEGEN) =>
           val toolConfig =
             if (options.contains('quantitative)) {
@@ -161,6 +169,8 @@ object KeYmaeraX {
       case Nil => map
       case "-help" :: _ => println(usage); exit(1)
       // actions
+      case "-taylorize" :: value :: tail =>
+        nextOption(map ++ Map('mode -> Modes.TAYLORIZE, 'file -> value), tail)
       case "-sandbox" :: tail =>
         nextOption(map ++ Map('sandbox -> true), tail)
       case "-modelplex" :: value :: tail =>
diff --git a/keymaerax-webui/src/main/scala/edu/cmu/cs/ls/keymaerax/launcher/TaylorizeMain.scala b/keymaerax-webui/src/main/scala/edu/cmu/cs/ls/keymaerax/launcher/TaylorizeMain.scala
new file mode 100644
index 0000000000..7d29c6a6cd
--- /dev/null
+++ b/keymaerax-webui/src/main/scala/edu/cmu/cs/ls/keymaerax/launcher/TaylorizeMain.scala
@@ -0,0 +1,63 @@
+/*
+ * Copyright (c) 2021 MIT-IBM Watson AI Lab, IBM Research.
+ */
+
+package edu.cmu.cs.ls.keymaerax.launcher
+
+import edu.cmu.cs.ls.keymaerax.btactics.helpers.DifferentialHelper
+import edu.cmu.cs.ls.keymaerax.core.{DifferentialProgram, ODESystem, PrettyPrinter, Term, Variable}
+import edu.cmu.cs.ls.keymaerax.parser.KeYmaeraXPrettyPrinter
+import edu.cmu.cs.ls.keymaerax.parser.StringConverter._
+
+
+/**
+  * The main class of Taylorizing sutff from the outside.
+  *
+  * @autor Nathan Fulton
+  */
+object TaylorizeMain {
+  PrettyPrinter.setPrinter(KeYmaeraXPrettyPrinter.pp)
+
+  def apply(fileName: String) = parseFile(fileName) match {
+    case Some(dp) => taylorize(dp) match {
+      case Some(listOfBounds) => listOfBounds.map(bounds => output(bounds._1, bounds._2, bounds._3)).reduce(_ + _)
+      case None => s"FAILED.\nKeYmaera X does not know how to construct a Taylor approximation of the system ${dp.prettyString}"
+    }
+    case None => s"FAILED.\nKeYmaera X does not know how to parse file ${fileName} with contents ${scala.io.Source.fromFile(fileName).mkString} into a DifferentialProgram."
+  }
+
+  def parseFile(fileName: String): Option[DifferentialProgram] = {
+    try {
+      val fileContents = scala.io.Source.fromFile(fileName).mkString
+      try {
+        Some(fileContents.asDifferentialProgram)
+      } catch {
+        case pe: edu.cmu.cs.ls.keymaerax.parser.ParseException => try {
+          Some(fileContents.asProgram.asInstanceOf[ODESystem].ode)
+        } catch {
+          case _: Throwable => None
+        }
+        case _: Throwable => None
+      }
+    }
+    catch {
+      case e: java.io.FileNotFoundException => None
+    }
+  }
+
+  def taylorize(dp: DifferentialProgram): Option[List[(Variable, Term, Term)]] = {
+    val vars = DifferentialHelper.atomicOdes(dp).map(atomic => atomic.xp.x)
+
+    //@todo Grab the taylor approximations from KeYmaera X backend FOR EACH PRIMED VARIABLE.
+    //@todo This whole section is complete nonsesnse filler.
+
+    Some(vars.map((v: Variable) => {
+      val name = v.name
+      (v, s"t^2/2".asTerm, s"t^3/3".asTerm)
+    }))
+  }
+
+  def output(v: Variable, lowerBound: Term, upperBound: Term) = {
+    s"\nvariablename:${v}\nlowerbound:${lowerBound.prettyString}\nupperbound:${upperBound.prettyString}\n"
+  }
+}

From cb18c9028bd45cce4d2060b38db7d9e6c1cd1dda Mon Sep 17 00:00:00 2001
From: Nathan Fulton <nathan@ibm.com>
Date: Wed, 30 Jun 2021 09:08:05 -0400
Subject: [PATCH 02/14] Improved error message in MathematicaLink.

The executor may be null if a tool or subtool isn't initialized before use.
This can happen, e.g., when adding a new subtool to the main Mathematica tool and
forgetting about the various pieces of boilerplate (calling .init, .shutdown, etc. in the
relevant methods of the Mathematica tool).

I've run into this issue a few times now, and I always lose a significant amount of time before
remember that tools have state and associated protocols. Hopefully this error message will make
tool extension/development easier.
---
 .../cmu/cs/ls/keymaerax/tools/ext/MathematicaLink.scala   | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/ext/MathematicaLink.scala b/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/ext/MathematicaLink.scala
index a219ec8600..868d2a1f46 100644
--- a/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/ext/MathematicaLink.scala
+++ b/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/ext/MathematicaLink.scala
@@ -313,6 +313,14 @@ class JLinkMathematicaLink(val engineName: String) extends MathematicaLink with
     */
   override def run[T](cmd: () => T, executor: ToolExecutor): T = {
     if (ml == null) throw new IllegalStateException("No MathKernel set")
+    if (executor == null) throw new IllegalStateException(
+      """
+        |No Executor was set.
+        |
+        |Likely explanation: a tool was used before initialization. Remember to call .init() on all tools or subtools
+        |before use. E.g., in edu.cmu.cs.ls.keymaerax.tools.ext.Mathematica, a call to .init() should be added for
+        |every new subtool.
+      """.stripMargin)
     val taskId = executor.schedule(_ => { ml.synchronized { cmd() } })
 
     executor.wait(taskId) match {

From da1300da3f40ee1c61561d6eb8db3fe892264c5f Mon Sep 17 00:00:00 2001
From: Nathan Fulton <nathan@ibm.com>
Date: Wed, 30 Jun 2021 09:11:45 -0400
Subject: [PATCH 03/14] Adds conversion support for Mathematica's C[].

Mathematica uses C[i] as its default form for parameters/constants.
For example, when computing series approxiamtions to solutions of ODEs, C[i]
is used to represent the constant of integration.

Perhaps a better way of handling this is to ensure that we never get C[i] back from
Mathematica. This can be achieved by always specifying initial conditions in calls to solvers.
For example, instead of:

```
AsymptoticDSolveValue[
  {
    x'[t] == y[t],
    y'[t] == x[t],
    y[0] == x0,
    x[0] == y0
  },
  {x[t], y[t]},
  {t, 0, 10}
]
```

instead do:

```
AsymptoticDSolveValue[
  {
    x'[t] == y[t],
    y'[t] == x[t],
    y[0] == x0,
    x[0] == y0
  },
  {x[t], y[t]},
  {t, 0, 10}
]
```

and then handle x0 and y0 before passing formulas back out of the tool.

Note: this is not currently used in any soundness-critical way that I'm aware of,
but both Reduce and DSolve are capable of generating C[i] symbols. Therefore,
we should do a thorough check that C[i]s are never handed back to us when not expected
or else revert this commit.

See https://reference.wolfram.com/language/ref/C.html for more information on what C[i]
is used for in Mathematica.
---
 .../cmu/cs/ls/keymaerax/tools/qe/MathematicaOpSpec.scala    | 2 ++
 .../cs/ls/keymaerax/tools/qe/MathematicaToKeYmaera.scala    | 6 ++++++
 2 files changed, 8 insertions(+)

diff --git a/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/qe/MathematicaOpSpec.scala b/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/qe/MathematicaOpSpec.scala
index 10f346d34c..7e5f84a582 100644
--- a/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/qe/MathematicaOpSpec.scala
+++ b/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/qe/MathematicaOpSpec.scala
@@ -129,6 +129,8 @@ object MathematicaOpSpec {
 
   def power: BinaryMathOpSpec = BinaryMathOpSpec(symbol("Power"))
 
+  def C: UnaryMathOpSpec = UnaryMathOpSpec(symbol("C")) //@todo document the meaning of this symbol.
+
   // implicit function application name[args]
   def func: NameMathOpSpec = NameMathOpSpec(
     (name: NamedSymbol, args: Array[Expr]) => {
diff --git a/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/qe/MathematicaToKeYmaera.scala b/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/qe/MathematicaToKeYmaera.scala
index cc2a5eb266..eec30bb6b0 100644
--- a/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/qe/MathematicaToKeYmaera.scala
+++ b/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/qe/MathematicaToKeYmaera.scala
@@ -41,6 +41,12 @@ class MathematicaToKeYmaera extends M2KConverter[KExpr] {
     //@note self-created MExpr with head RATIONAL are not rationalQ (type identifiers do not match)
     else if (MathematicaOpSpec.rational.applies(e)) convertBinary(e, Divide.apply)
 
+    // Constant symbols, typically as constants of integration.
+    // Should NOT be leaked out of tooling code. Marking as "interpreted" because that should prevent its unsound use in
+    // at least the most important soundness-critical contexts.
+    // @todo either enforce this invariant or give the function an obviously obnoxious name as a warning to the user.
+    else if (MathematicaOpSpec.C.applies(e))     convertUnary (e, (t: Term) => FuncOf(Function("C", None, Real, Real, true), t))
+
     // Arith expressions
     else if (MathematicaOpSpec.plus.applies(e))   convertNary  (e, Plus.apply)
     else if (MathematicaOpSpec.minus.applies(e))  convertBinary(e, Minus.apply)

From db82009d5a97bc181d1a3261a4f89c79b6dfab55 Mon Sep 17 00:00:00 2001
From: Nathan Fulton <nathan@ibm.com>
Date: Wed, 30 Jun 2021 09:18:12 -0400
Subject: [PATCH 04/14] Convert DifferentialPrograms to Mathematica.

Note: this code would already written in the DSolve tool, and now
exists in two different places. We should unify this code and clean up
the interface to the DifferenitalProgram converter.
---
 .../qe/DifferentialProgramToMathematica.scala | 81 +++++++++++++++++++
 1 file changed, 81 insertions(+)
 create mode 100644 keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/qe/DifferentialProgramToMathematica.scala

diff --git a/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/qe/DifferentialProgramToMathematica.scala b/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/qe/DifferentialProgramToMathematica.scala
new file mode 100644
index 0000000000..89a0441974
--- /dev/null
+++ b/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/qe/DifferentialProgramToMathematica.scala
@@ -0,0 +1,81 @@
+package edu.cmu.cs.ls.keymaerax.tools.qe
+
+import edu.cmu.cs.ls.keymaerax.core.{AtomicODE, BaseVariable, DifferentialProduct, DifferentialProgram, DifferentialSymbol, Equal, FuncOf, Function, Number, ODESystem, Program, Term, Variable}
+import edu.cmu.cs.ls.keymaerax.infrastruct.ExpressionTraversal.{ExpressionTraversalFunction, StopTraversal}
+import edu.cmu.cs.ls.keymaerax.infrastruct.{ExpressionTraversal, PosInExpr}
+import edu.cmu.cs.ls.keymaerax.tools.ConversionException
+import edu.cmu.cs.ls.keymaerax.tools.ext.ExtMathematicaOpSpec
+import edu.cmu.cs.ls.keymaerax.tools.qe.MathematicaConversion.MExpr
+
+import scala.collection.immutable.{List, Map}
+import scala.math.BigDecimal
+
+case class DifferentialProgramToMathematica(k2m: K2MConverter[MathematicaConversion.KExpr]) {
+  /**
+    * Converts a DifferentialProgram to a Mathematica expression.
+    * This code is largely copy/pasted from MathematicaTools.MathematicalODESolvertool, so that we can KeYmaera
+    * translations of ODEs in other contexts.
+    * @return a 4-tuple containing:
+    *         a list of MExprs specifying the ODEs,
+    *         a list of MExprs specifying symbolic initial conditions,
+    *         a list of MExprs specifying functions (IDK what this is actually -- check the dsolve code and Wolfram documentaiton to figure out)
+    *         the variable used as the time variable (diffArg, as passed in initially)
+    * @author Nathan Fulton
+    */
+  def apply(dp: DifferentialProgram, diffArg: Variable, iv: Map[Variable, Variable]): (List[MExpr], List[MExpr], List[MExpr], Variable) = {
+    /**
+      * @note copied from MathematicaTools.MathematicalODESolvertool.
+      */
+    def toDiffSys(diffSys: DifferentialProgram, diffArg: Variable): List[(Variable, Term)] = {
+      var result = List[(Variable, Term)]()
+      ExpressionTraversal.traverse(new ExpressionTraversalFunction {
+        override def preP(p: PosInExpr, e: Program): Either[Option[StopTraversal], Program] = e match {
+          case AtomicODE(DifferentialSymbol(x), theta) if x != diffArg => result = result :+ (x, theta); Left(None)
+          case AtomicODE(DifferentialSymbol(x), _) if x == diffArg => Left(None)
+          case ODESystem(_, _) => Left(None)
+          case DifferentialProduct(_, _) => Left(None)
+        }
+      }, diffSys)
+      result
+    }
+
+    /** @note coped from dsolve tool. */
+    def functionalizeVars(t: Term, arg: Term, vars: Variable*) = ExpressionTraversal.traverse(
+      new ExpressionTraversalFunction {
+        override def postT(p: PosInExpr, e: Term): Either[Option[StopTraversal], Term] = e match {
+          case v@BaseVariable(name, idx, sort) if vars.isEmpty || vars.contains(v) =>
+            Right(FuncOf(Function(name, idx, arg.sort, sort), arg))
+          case _ => Left(None)
+        }
+      }, t) match {
+      case Some(resultTerm) => resultTerm
+      case None => throw ConversionException("Unable to functionalize " + t)
+    }
+
+    val diffSys = toDiffSys(dp, diffArg)
+
+    val primedVars = diffSys.map(_._1)
+    val functionalizedTerms = diffSys.map{ case (x, theta) => ( x, functionalizeVars(theta, diffArg, primedVars:_*)) }
+    val mathTerms = functionalizedTerms.map({case (x, theta) =>
+      (ExtMathematicaOpSpec.dx(ExtMathematicaOpSpec.primed(k2m(x)))(k2m(diffArg)), k2m(theta))})
+    val convertedDiffSys = mathTerms.map({case (x, theta) => MathematicaOpSpec.equal(x, theta)})
+
+    val functions = diffSys.map(t => k2m(functionalizeVars(t._1, diffArg)))
+
+    //@todo allows for partial initial conditions, but this will probably cause issues if IVs are used.
+    val initialValues = diffSys
+      .map(t =>
+        iv.get(t._1) match {
+          case Some(definedInitialValue) =>
+            Some(
+              Equal(functionalizeVars(t._1, Number(BigDecimal(0)), primedVars:_*), definedInitialValue)
+            )
+          case None => None //@todo perhaps add an error here?
+        }
+      )
+      .filterNot(x => x.isEmpty)
+      .map(x => k2m(x.get))
+
+    (convertedDiffSys, initialValues, functions, diffArg)
+  }
+}

From e5cc95a9a4e579a18e164f45c11b428e3b3a318e Mon Sep 17 00:00:00 2001
From: Nathan Fulton <nathan@ibm.com>
Date: Wed, 30 Jun 2021 09:26:18 -0400
Subject: [PATCH 05/14] Tool support for Mathematica's AsymptoticDSolveValue

Computes asymptotic approximations to the solutions of ODEs.

See https://reference.wolfram.com/language/ref/AsymptoticDSolveValue.html
---
 ...ntialSolutionSeriesApproximationTool.scala | 14 ++++
 .../tools/ext/ExtMathematicaOpSpec.scala      |  2 +
 .../ls/keymaerax/tools/ext/Mathematica.scala  | 12 ++-
 .../tools/ext/MathematicaTools.scala          | 75 ++++++++++++++++++-
 .../ls/keymaerax/launcher/TaylorizeMain.scala |  4 +-
 .../tools/SeriesExpansionTests.scala          | 31 ++++++++
 6 files changed, 133 insertions(+), 5 deletions(-)
 create mode 100644 keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/ext/DifferentialSolutionSeriesApproximationTool.scala
 create mode 100644 keymaerax-webui/src/test/scala/edu/cmu/cs/ls/keymaerax/tools/SeriesExpansionTests.scala

diff --git a/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/ext/DifferentialSolutionSeriesApproximationTool.scala b/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/ext/DifferentialSolutionSeriesApproximationTool.scala
new file mode 100644
index 0000000000..0d6697324f
--- /dev/null
+++ b/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/ext/DifferentialSolutionSeriesApproximationTool.scala
@@ -0,0 +1,14 @@
+package edu.cmu.cs.ls.keymaerax.tools.ext
+
+import edu.cmu.cs.ls.keymaerax.core.{NamedSymbol, Number, ODESystem, Term, Variable}
+import edu.cmu.cs.ls.keymaerax.tools.ToolInterface
+
+trait DifferentialSolutionSeriesApproximationTool extends ToolInterface {
+  /**
+    *
+    * @param odes The ODEs
+    * @param ctx Context for any parameters in the ODEs.
+    * @return upper and lower bound series approximations of each primed variable in the ODEs.
+    */
+  def seriesApproximation(odes: ODESystem, ctx: Map[Term, Term]): Option[Map[Variable, Term]]
+}
diff --git a/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/ext/ExtMathematicaOpSpec.scala b/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/ext/ExtMathematicaOpSpec.scala
index 4c654734eb..9e02249696 100644
--- a/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/ext/ExtMathematicaOpSpec.scala
+++ b/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/ext/ExtMathematicaOpSpec.scala
@@ -101,6 +101,8 @@ object ExtMathematicaOpSpec {
 
   def d: BinaryMathOpSpec = BinaryMathOpSpec(symbol("D"))
 
+  def dsolveAsymptoticApproximation: NaryMathOpSpec = NaryMathOpSpec(symbol("AsymptoticDSolveValue"))
+
   //</editor-fold>
 
 }
diff --git a/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/ext/Mathematica.scala b/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/ext/Mathematica.scala
index f89b131ef9..187ea34bcd 100644
--- a/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/ext/Mathematica.scala
+++ b/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/ext/Mathematica.scala
@@ -15,6 +15,7 @@ import edu.cmu.cs.ls.keymaerax.pt.ProvableSig
 import edu.cmu.cs.ls.keymaerax.tools.ext.SimulationTool.{SimRun, SimState, Simulation}
 import edu.cmu.cs.ls.keymaerax.tools._
 import edu.cmu.cs.ls.keymaerax.tools.ext.SOSsolveTool.Result
+import edu.cmu.cs.ls.keymaerax.tools.ext.MathematicaDifferentialSolutionSeriesApproximationTool
 
 import scala.annotation.tailrec
 import scala.collection.immutable.{Map, Seq}
@@ -30,7 +31,7 @@ import scala.collection.immutable.{Map, Seq}
 class Mathematica(private[tools] val link: MathematicaLink, override val name: String) extends Tool
     with QETacticTool with InvGenTool with ODESolverTool with CounterExampleTool
     with SimulationTool with DerivativeTool with EquationSolverTool with SimplificationTool with AlgebraTool
-    with PDESolverTool with SOSsolveTool with ToolOperationManagement {
+    with PDESolverTool with SOSsolveTool with ToolOperationManagement with DifferentialSolutionSeriesApproximationTool {
 
   /** Indicates whether the tool is initialized. */
   private var initialized = false
@@ -45,6 +46,7 @@ class Mathematica(private[tools] val link: MathematicaLink, override val name: S
   private val mSolve = new MathematicaEquationSolverTool(link)
   private val mAlgebra = new MathematicaAlgebraTool(link)
   private val mSimplify = new MathematicaSimplificationTool(link)
+  private val mExpand = new MathematicaDifferentialSolutionSeriesApproximationTool(link)
 
   private val qeInitialTimeout = Integer.parseInt(Configuration(Configuration.Keys.QE_TIMEOUT_INITIAL))
   private val qeCexTimeout = Integer.parseInt(Configuration(Configuration.Keys.QE_TIMEOUT_CEX))
@@ -63,6 +65,7 @@ class Mathematica(private[tools] val link: MathematicaLink, override val name: S
     mSolve.memoryLimit = memoryLimit
     mAlgebra.memoryLimit = memoryLimit
     mSimplify.memoryLimit = memoryLimit
+    mExpand.memoryLimit = memoryLimit
 
     // initialze tool thread pools
     mQE.init()
@@ -75,6 +78,8 @@ class Mathematica(private[tools] val link: MathematicaLink, override val name: S
     mAlgebra.init()
     mSimplify.init()
     mSOSsolve.init()
+    mExpand.init()
+
 
     initialized = link match {
       case l: JLinkMathematicaLink =>
@@ -106,6 +111,7 @@ class Mathematica(private[tools] val link: MathematicaLink, override val name: S
     mSolve.shutdown()
     mAlgebra.shutdown()
     mSimplify.shutdown()
+    mExpand.shutdown()
     //@note last, because we want to shut down all executors (tool threads) before shutting down the JLink interface
     link match {
       case l: JLinkMathematicaLink => l.shutdown()
@@ -230,4 +236,8 @@ class Mathematica(private[tools] val link: MathematicaLink, override val name: S
 
   /** @inheritdoc */
   override def getAvailableWorkers: Int = 1
+
+  /** @inheritdoc */
+  override def seriesApproximation(odes: ODESystem, ctx: Map[Term, Term]): Option[Map[Variable, Term]] =
+    mExpand.seriesApproximation(odes, ctx)
 }
diff --git a/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/ext/MathematicaTools.scala b/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/ext/MathematicaTools.scala
index 44d9ac048a..54d4b1592b 100644
--- a/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/ext/MathematicaTools.scala
+++ b/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/ext/MathematicaTools.scala
@@ -10,12 +10,13 @@ import edu.cmu.cs.ls.keymaerax.btactics.InvariantGenerator
 import edu.cmu.cs.ls.keymaerax.btactics.helpers.DifferentialHelper
 import edu.cmu.cs.ls.keymaerax.core.{Variable, _}
 import edu.cmu.cs.ls.keymaerax.infrastruct.ExpressionTraversal.{ExpressionTraversalFunction, StopTraversal}
-import edu.cmu.cs.ls.keymaerax.infrastruct.{ExpressionTraversal, FormulaTools, PosInExpr}
+import edu.cmu.cs.ls.keymaerax.infrastruct.{ExpressionTraversal, FormulaTools, PosInExpr, SubstitutionHelper}
 import edu.cmu.cs.ls.keymaerax.tools.qe.MathematicaConversion.{KExpr, _}
 import edu.cmu.cs.ls.keymaerax.tools.ext.SimulationTool.{SimRun, SimState, Simulation}
-import edu.cmu.cs.ls.keymaerax.tools.qe.{BinaryMathOpSpec, ExprFactory, K2MConverter, KeYmaeraToMathematica, M2KConverter, MathematicaNameConversion, MathematicaOpSpec, MathematicaToKeYmaera, NaryMathOpSpec, UnaryMathOpSpec}
+import edu.cmu.cs.ls.keymaerax.tools.qe.{DifferentialProgramToMathematica, BinaryMathOpSpec, ExprFactory, K2MConverter, KeYmaeraToMathematica, M2KConverter, MathematicaNameConversion, MathematicaOpSpec, MathematicaToKeYmaera, NaryMathOpSpec, UnaryMathOpSpec}
 import edu.cmu.cs.ls.keymaerax.tools.qe.MathematicaOpSpec._
 import edu.cmu.cs.ls.keymaerax.tools._
+import edu.cmu.cs.ls.keymaerax.infrastruct.Augmentors._
 
 import scala.collection.immutable
 import scala.math.BigDecimal
@@ -755,6 +756,76 @@ class MathematicaEquationSolverTool(override val link: MathematicaLink) extends
   }
 }
 
+/**
+  * Uses Mathematica's AsymptoticDSolveValue function to construct power series expansions of the solution to a system of ODEs.
+  *
+  * @author Nathan Fulton
+  */
+class MathematicaDifferentialSolutionSeriesApproximationTool(override val link: MathematicaLink) extends BaseKeYmaeraMathematicaBridge[KExpr](link, new UncheckedBaseK2MConverter, PegasusM2KConverter) with DifferentialSolutionSeriesApproximationTool {
+
+  private def pairToList(t: Term): List[Term] = t match {
+    case Pair(left, right) => {
+      assert(!left.isInstanceOf[Pair], "Expected canonical ordering so that tail recursion is possible.")
+      left :: pairToList(right)
+    }
+    case _ => t :: Nil
+  }
+
+  private def convertResult(odes: ODESystem, result: Expression) = {
+    result match {
+      case t: Term => {
+        // Convert nested pairs into a list of terms.
+        val approximations = pairToList(t)
+
+        // Replace all constant symbols in the terms with old(x_i) where x_i is the i^th primed variable in the ODEs.
+        val primedVariables = DifferentialHelper.getPrimedVariables(odes)
+        val approximationsWithInitialConditions = approximations.map(approximation => {
+          Range(1, primedVariables.length + 1).foldLeft(approximation)((approximation, i) => {
+            val constantSymbol = FuncOf(Function("C", None, Real, Real, true), Number(i)) // C[i]
+            val x_i = primedVariables(i - 1) //the i^th primed variable in the system.
+            val old_x_i = FuncOf(Function("old", None, Real, Real), x_i) // old(x_i)
+            approximation.replaceAll(constantSymbol, old_x_i)
+          })
+        })
+
+        Some(
+          DifferentialHelper.getPrimedVariables(odes)
+            .zip(approximationsWithInitialConditions)
+            .toMap
+        )
+      }
+      case _ => throw new ConversionException(s"Expected AsymptoticDSolveValue to give a term or a list of terms, but found ${result}")
+    }
+  }
+
+  /** @inheritdoc */
+  override def seriesApproximation(odes: ODESystem, ctx: Map[Term, Term]): Option[Map[Variable, Term]] = {
+    //See https://reference.wolfram.com/language/ref/AsymptoticDSolveValue.html for the command we're going to use.
+
+    // Apply the context to the ODEs.
+    val odesInCtx = ctx.foldLeft(odes)((currExpr, x) => {
+      SubstitutionHelper.replaceFree(currExpr)(x._1, x._2).asInstanceOf[ODESystem]
+    })
+
+    // Identify a unique time variable.
+    val TIME_VAR = Variable("t") //@todo do something here that's less stupid.
+
+    // Convert the ODEs into Mathematica expressions.
+    // @todo move this out into a utility class independent of the k2m converter interface.
+    val (mODEs, ivs, _, _) = DifferentialProgramToMathematica(k2m).apply(odesInCtx.ode, TIME_VAR, Map())
+
+    //construct the Mathematica expression that will execute AsymptoticDSolveValue.
+    val input = ExtMathematicaOpSpec.dsolveAsymptoticApproximation(
+      MathematicaOpSpec.list(mODEs:_*),
+      MathematicaOpSpec.list(DifferentialHelper.getPrimedVariables(odesInCtx).map(k2m):_*),
+      MathematicaOpSpec.list(k2m(TIME_VAR), k2m(Number(0)), k2m(Number(10))) //@todo I have no clue what these bounds should be. Make these arguments at least.
+    )
+
+    val (_, result) = run(input)
+    convertResult(odes, result)
+  }
+}
+
 /**
   * A link to Mathematica using the JLink interface.
   *
diff --git a/keymaerax-webui/src/main/scala/edu/cmu/cs/ls/keymaerax/launcher/TaylorizeMain.scala b/keymaerax-webui/src/main/scala/edu/cmu/cs/ls/keymaerax/launcher/TaylorizeMain.scala
index 7d29c6a6cd..c542af2436 100644
--- a/keymaerax-webui/src/main/scala/edu/cmu/cs/ls/keymaerax/launcher/TaylorizeMain.scala
+++ b/keymaerax-webui/src/main/scala/edu/cmu/cs/ls/keymaerax/launcher/TaylorizeMain.scala
@@ -11,9 +11,9 @@ import edu.cmu.cs.ls.keymaerax.parser.StringConverter._
 
 
 /**
-  * The main class of Taylorizing sutff from the outside.
+  * The main class for the CLI interface to Taylor approximating solutions to ODEs.
   *
-  * @autor Nathan Fulton
+  * @author Nathan Fulton
   */
 object TaylorizeMain {
   PrettyPrinter.setPrinter(KeYmaeraXPrettyPrinter.pp)
diff --git a/keymaerax-webui/src/test/scala/edu/cmu/cs/ls/keymaerax/tools/SeriesExpansionTests.scala b/keymaerax-webui/src/test/scala/edu/cmu/cs/ls/keymaerax/tools/SeriesExpansionTests.scala
new file mode 100644
index 0000000000..462775ad68
--- /dev/null
+++ b/keymaerax-webui/src/test/scala/edu/cmu/cs/ls/keymaerax/tools/SeriesExpansionTests.scala
@@ -0,0 +1,31 @@
+package edu.cmu.cs.ls.keymaerax.tools
+
+import edu.cmu.cs.ls.keymaerax.btactics.TacticTestBase
+import edu.cmu.cs.ls.keymaerax.core.{NamedSymbol, Number, ODESystem, Term}
+import edu.cmu.cs.ls.keymaerax.parser.StringConverter._
+import smtlib.theories.Core.True
+
+class SeriesExpansionTests extends TacticTestBase {
+  "series expansion tool" should "work for a simple system" in withMathematica(tool => {
+    val odes = ODESystem("{x'=y, y'=-a*x}".asDifferentialProgram, "true".asFormula)
+    val ctx: Map[Term, Term] = Map({"a".asVariable -> Number(5)})
+    val result = tool.seriesApproximation(odes, ctx)
+    println(result)
+  })
+
+  it should "return some result when there's no relevant context." in withMathematica(tool => {
+    val odes = ODESystem("{x'=y, y'=-x}".asDifferentialProgram, "true".asFormula)
+    val ctx: Map[Term, Term] = Map()
+    val result = tool.seriesApproximation(odes, ctx)
+    println(result)
+  })
+
+  it should "return the correct result for x'=1" in withMathematica(tool => {
+    val odes = ODESystem("x'=1".asDifferentialProgram, "true".asFormula)
+    val ctx: Map[Term, Term] = Map()
+    val result = tool.seriesApproximation(odes, ctx)
+
+    result.get("x".asVariable) shouldBe "t + old(x)".asTerm
+  })
+
+}

From 0e88b6a2a713281ace19462edfce0252a7b3a83c Mon Sep 17 00:00:00 2001
From: Nathan Fulton <nathan@ibm.com>
Date: Wed, 30 Jun 2021 10:00:26 -0400
Subject: [PATCH 06/14] adds a default tool provider for series approximations.

---
 .../edu/cmu/cs/ls/keymaerax/btactics/ToolProvider.scala   | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/btactics/ToolProvider.scala b/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/btactics/ToolProvider.scala
index 4e86c2255d..0999007545 100644
--- a/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/btactics/ToolProvider.scala
+++ b/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/btactics/ToolProvider.scala
@@ -68,6 +68,8 @@ object ToolProvider extends ToolProvider with Logging {
 
   def sosSolveTool(): Option[SOSsolveTool] = f.sosSolveTool()
 
+  def differentialSeriesApproxmationnTool() = f.differentialSeriesApproxmationnTool()
+
   def init(): Boolean = f.init()
 
   def shutdown(): Unit = f.shutdown()
@@ -147,6 +149,9 @@ trait ToolProvider {
   /** Returns a SOSsolve tool. */
   def sosSolveTool(): Option[SOSsolveTool]
 
+  /** Returns a series expansion tool. */
+  def differentialSeriesApproxmationnTool(): Option[DifferentialSolutionSeriesApproximationTool]
+
   /** Initializes the tools. */
   def init(): Boolean
 
@@ -174,6 +179,7 @@ class PreferredToolProvider[T <: Tool](val toolPreferences: List[T]) extends Too
   private[this] lazy val solver: Option[Tool with EquationSolverTool] = toolPreferences.find(_.isInstanceOf[EquationSolverTool]).map(_.asInstanceOf[Tool with EquationSolverTool])
   private[this] lazy val algebra: Option[Tool with AlgebraTool] = toolPreferences.find(_.isInstanceOf[AlgebraTool]).map(_.asInstanceOf[Tool with AlgebraTool])
   private[this] lazy val sossolve: Option[Tool with SOSsolveTool] = toolPreferences.find(_.isInstanceOf[SOSsolveTool]).map(_.asInstanceOf[Tool with SOSsolveTool])
+  private[this] lazy val diffSeriesApproximation: Option[Tool with DifferentialSolutionSeriesApproximationTool] = toolPreferences.find(_.isInstanceOf[DifferentialSolutionSeriesApproximationTool]).map(_.asInstanceOf[Tool with DifferentialSolutionSeriesApproximationTool])
 
   override def tools(): List[Tool] = toolPreferences
   override def defaultTool(): Option[Tool] = toolPreferences.headOption
@@ -193,6 +199,7 @@ class PreferredToolProvider[T <: Tool](val toolPreferences: List[T]) extends Too
   override def solverTool(): Option[EquationSolverTool] = ensureInitialized(solver)
   override def algebraTool(): Option[AlgebraTool] = ensureInitialized(algebra)
   override def sosSolveTool(): Option[SOSsolveTool] = ensureInitialized(sossolve)
+  override def differentialSeriesApproxmationnTool(): Option[DifferentialSolutionSeriesApproximationTool] = ensureInitialized(diffSeriesApproximation)
   override def init(): Boolean = false /* override to initialize tools in more specialized providers */
   override def shutdown(): Unit = toolPreferences.foreach(_.shutdown())
   override def isInitialized: Boolean = toolPreferences.forall(_.isInitialized)
@@ -221,6 +228,7 @@ class NoneToolProvider extends ToolProvider {
   override def solverTool(): Option[EquationSolverTool] = None
   override def algebraTool(): Option[AlgebraTool] = None
   override def sosSolveTool(): Option[SOSsolveTool] = None
+  override def differentialSeriesApproxmationnTool() = None
   override def init(): Boolean = true
   override def shutdown(): Unit = {}
   override def isInitialized: Boolean = true

From 4722aa2a62f7dac32d90baab66847f49dbd4f705 Mon Sep 17 00:00:00 2001
From: Nathan Fulton <nathan@ibm.com>
Date: Wed, 30 Jun 2021 10:01:01 -0400
Subject: [PATCH 07/14] incoroproate new series approximation tool into the
 (perhaps poorly named) taylorize cli mode.

---
 .../cmu/cs/ls/keymaerax/launcher/KeYmaeraX.scala |  8 ++++++++
 .../cs/ls/keymaerax/launcher/TaylorizeMain.scala | 16 ++++++++++------
 2 files changed, 18 insertions(+), 6 deletions(-)

diff --git a/keymaerax-webui/src/main/scala/edu/cmu/cs/ls/keymaerax/launcher/KeYmaeraX.scala b/keymaerax-webui/src/main/scala/edu/cmu/cs/ls/keymaerax/launcher/KeYmaeraX.scala
index d8f7b60a3a..bf29c75485 100644
--- a/keymaerax-webui/src/main/scala/edu/cmu/cs/ls/keymaerax/launcher/KeYmaeraX.scala
+++ b/keymaerax-webui/src/main/scala/edu/cmu/cs/ls/keymaerax/launcher/KeYmaeraX.scala
@@ -97,6 +97,14 @@ object KeYmaeraX {
       //@todo allow multiple passes by filter architecture: -prove bla.key -tactic bla.scal -modelplex -codegen
       options.get('mode) match {
         case Some(Modes.TAYLORIZE) => {
+          val toolConfig =
+            if (options.contains('quantitative)) {
+              configFromFile(Tools.MATHEMATICA) //@note quantitative ModelPlex uses Mathematica to simplify formulas
+            } else {
+              configFromFile("z3")
+            }
+          initializeProver(combineConfigs(options, toolConfig), usage)
+
           val filename = options.get('file)
           filename match {
             case Some(s) => println(TaylorizeMain(s.asInstanceOf[String]))
diff --git a/keymaerax-webui/src/main/scala/edu/cmu/cs/ls/keymaerax/launcher/TaylorizeMain.scala b/keymaerax-webui/src/main/scala/edu/cmu/cs/ls/keymaerax/launcher/TaylorizeMain.scala
index c542af2436..beef72d66f 100644
--- a/keymaerax-webui/src/main/scala/edu/cmu/cs/ls/keymaerax/launcher/TaylorizeMain.scala
+++ b/keymaerax-webui/src/main/scala/edu/cmu/cs/ls/keymaerax/launcher/TaylorizeMain.scala
@@ -4,6 +4,8 @@
 
 package edu.cmu.cs.ls.keymaerax.launcher
 
+import edu.cmu.cs.ls.keymaerax.btactics.ToolProvider
+import edu.cmu.cs.ls.keymaerax.{Configuration, FileConfiguration}
 import edu.cmu.cs.ls.keymaerax.btactics.helpers.DifferentialHelper
 import edu.cmu.cs.ls.keymaerax.core.{DifferentialProgram, ODESystem, PrettyPrinter, Term, Variable}
 import edu.cmu.cs.ls.keymaerax.parser.KeYmaeraXPrettyPrinter
@@ -48,13 +50,15 @@ object TaylorizeMain {
   def taylorize(dp: DifferentialProgram): Option[List[(Variable, Term, Term)]] = {
     val vars = DifferentialHelper.atomicOdes(dp).map(atomic => atomic.xp.x)
 
-    //@todo Grab the taylor approximations from KeYmaera X backend FOR EACH PRIMED VARIABLE.
-    //@todo This whole section is complete nonsesnse filler.
+    val result = ToolProvider.differentialSeriesApproxmationnTool().get.seriesApproximation(
+      ODESystem(dp, "true".asFormula),
+      Map()
+    )
 
-    Some(vars.map((v: Variable) => {
-      val name = v.name
-      (v, s"t^2/2".asTerm, s"t^3/3".asTerm)
-    }))
+    result match {
+      case Some(mapping) => Some(mapping.map(vt => (vt._1,vt._2,vt._2)).toList) //@todo get upper and lower bounds.
+      case None => None
+    }
   }
 
   def output(v: Variable, lowerBound: Term, upperBound: Term) = {

From 7d73eab79b634287f5d74199bfab5fa4be08df0c Mon Sep 17 00:00:00 2001
From: Nathan Fulton <nathan@ibm.com>
Date: Tue, 22 Jun 2021 12:08:44 -0400
Subject: [PATCH 08/14] taylorize interface

---
 .../cs/ls/keymaerax/launcher/KeYmaeraX.scala  | 12 +++-
 .../ls/keymaerax/launcher/TaylorizeMain.scala | 63 +++++++++++++++++++
 2 files changed, 74 insertions(+), 1 deletion(-)
 create mode 100644 keymaerax-webui/src/main/scala/edu/cmu/cs/ls/keymaerax/launcher/TaylorizeMain.scala

diff --git a/keymaerax-webui/src/main/scala/edu/cmu/cs/ls/keymaerax/launcher/KeYmaeraX.scala b/keymaerax-webui/src/main/scala/edu/cmu/cs/ls/keymaerax/launcher/KeYmaeraX.scala
index df3f018605..6d6a2f6fc6 100644
--- a/keymaerax-webui/src/main/scala/edu/cmu/cs/ls/keymaerax/launcher/KeYmaeraX.scala
+++ b/keymaerax-webui/src/main/scala/edu/cmu/cs/ls/keymaerax/launcher/KeYmaeraX.scala
@@ -67,7 +67,8 @@ object KeYmaeraX {
     val CONVERT: String = edu.cmu.cs.ls.keymaerax.cli.KeYmaeraX.Modes.CONVERT
     val SETUP: String = edu.cmu.cs.ls.keymaerax.cli.KeYmaeraX.Modes.SETUP
     val UI: String = "ui"
-    val modes: Set[String] = Set(CODEGEN, CONVERT, MODELPLEX, PROVE, REPL, UI, SETUP)
+    val TAYLORIZE: String = "taylorize"
+    val modes: Set[String] = Set(CODEGEN, CONVERT, MODELPLEX, PROVE, REPL, UI, SETUP, TAYLORIZE)
   }
 
   /** Usage -help information. */
@@ -95,6 +96,13 @@ object KeYmaeraX {
     try {
       //@todo allow multiple passes by filter architecture: -prove bla.key -tactic bla.scal -modelplex -codegen
       options.get('mode) match {
+        case Some(Modes.TAYLORIZE) => {
+          val filename = options.get('file)
+          filename match {
+            case Some(s) => println(TaylorizeMain(s.asInstanceOf[String]))
+            case None => println("FAILED.")
+          }
+        }
         case Some(Modes.CODEGEN) =>
           val toolConfig =
             if (options.contains('quantitative)) {
@@ -161,6 +169,8 @@ object KeYmaeraX {
       case Nil => map
       case "-help" :: _ => println(usage); exit(1)
       // actions
+      case "-taylorize" :: value :: tail =>
+        nextOption(map ++ Map('mode -> Modes.TAYLORIZE, 'file -> value), tail)
       case "-sandbox" :: tail =>
         nextOption(map ++ Map('sandbox -> true), tail)
       case "-modelplex" :: value :: tail =>
diff --git a/keymaerax-webui/src/main/scala/edu/cmu/cs/ls/keymaerax/launcher/TaylorizeMain.scala b/keymaerax-webui/src/main/scala/edu/cmu/cs/ls/keymaerax/launcher/TaylorizeMain.scala
new file mode 100644
index 0000000000..7d29c6a6cd
--- /dev/null
+++ b/keymaerax-webui/src/main/scala/edu/cmu/cs/ls/keymaerax/launcher/TaylorizeMain.scala
@@ -0,0 +1,63 @@
+/*
+ * Copyright (c) 2021 MIT-IBM Watson AI Lab, IBM Research.
+ */
+
+package edu.cmu.cs.ls.keymaerax.launcher
+
+import edu.cmu.cs.ls.keymaerax.btactics.helpers.DifferentialHelper
+import edu.cmu.cs.ls.keymaerax.core.{DifferentialProgram, ODESystem, PrettyPrinter, Term, Variable}
+import edu.cmu.cs.ls.keymaerax.parser.KeYmaeraXPrettyPrinter
+import edu.cmu.cs.ls.keymaerax.parser.StringConverter._
+
+
+/**
+  * The main class of Taylorizing sutff from the outside.
+  *
+  * @autor Nathan Fulton
+  */
+object TaylorizeMain {
+  PrettyPrinter.setPrinter(KeYmaeraXPrettyPrinter.pp)
+
+  def apply(fileName: String) = parseFile(fileName) match {
+    case Some(dp) => taylorize(dp) match {
+      case Some(listOfBounds) => listOfBounds.map(bounds => output(bounds._1, bounds._2, bounds._3)).reduce(_ + _)
+      case None => s"FAILED.\nKeYmaera X does not know how to construct a Taylor approximation of the system ${dp.prettyString}"
+    }
+    case None => s"FAILED.\nKeYmaera X does not know how to parse file ${fileName} with contents ${scala.io.Source.fromFile(fileName).mkString} into a DifferentialProgram."
+  }
+
+  def parseFile(fileName: String): Option[DifferentialProgram] = {
+    try {
+      val fileContents = scala.io.Source.fromFile(fileName).mkString
+      try {
+        Some(fileContents.asDifferentialProgram)
+      } catch {
+        case pe: edu.cmu.cs.ls.keymaerax.parser.ParseException => try {
+          Some(fileContents.asProgram.asInstanceOf[ODESystem].ode)
+        } catch {
+          case _: Throwable => None
+        }
+        case _: Throwable => None
+      }
+    }
+    catch {
+      case e: java.io.FileNotFoundException => None
+    }
+  }
+
+  def taylorize(dp: DifferentialProgram): Option[List[(Variable, Term, Term)]] = {
+    val vars = DifferentialHelper.atomicOdes(dp).map(atomic => atomic.xp.x)
+
+    //@todo Grab the taylor approximations from KeYmaera X backend FOR EACH PRIMED VARIABLE.
+    //@todo This whole section is complete nonsesnse filler.
+
+    Some(vars.map((v: Variable) => {
+      val name = v.name
+      (v, s"t^2/2".asTerm, s"t^3/3".asTerm)
+    }))
+  }
+
+  def output(v: Variable, lowerBound: Term, upperBound: Term) = {
+    s"\nvariablename:${v}\nlowerbound:${lowerBound.prettyString}\nupperbound:${upperBound.prettyString}\n"
+  }
+}

From 740837feca1027b86c6bbf3f2320e20cf6f4dcc1 Mon Sep 17 00:00:00 2001
From: Nathan Fulton <nathan@ibm.com>
Date: Wed, 30 Jun 2021 09:08:05 -0400
Subject: [PATCH 09/14] Improved error message in MathematicaLink.

The executor may be null if a tool or subtool isn't initialized before use.
This can happen, e.g., when adding a new subtool to the main Mathematica tool and
forgetting about the various pieces of boilerplate (calling .init, .shutdown, etc. in the
relevant methods of the Mathematica tool).

I've run into this issue a few times now, and I always lose a significant amount of time before
remember that tools have state and associated protocols. Hopefully this error message will make
tool extension/development easier.
---
 .../cmu/cs/ls/keymaerax/tools/ext/MathematicaLink.scala   | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/ext/MathematicaLink.scala b/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/ext/MathematicaLink.scala
index a219ec8600..868d2a1f46 100644
--- a/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/ext/MathematicaLink.scala
+++ b/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/ext/MathematicaLink.scala
@@ -313,6 +313,14 @@ class JLinkMathematicaLink(val engineName: String) extends MathematicaLink with
     */
   override def run[T](cmd: () => T, executor: ToolExecutor): T = {
     if (ml == null) throw new IllegalStateException("No MathKernel set")
+    if (executor == null) throw new IllegalStateException(
+      """
+        |No Executor was set.
+        |
+        |Likely explanation: a tool was used before initialization. Remember to call .init() on all tools or subtools
+        |before use. E.g., in edu.cmu.cs.ls.keymaerax.tools.ext.Mathematica, a call to .init() should be added for
+        |every new subtool.
+      """.stripMargin)
     val taskId = executor.schedule(_ => { ml.synchronized { cmd() } })
 
     executor.wait(taskId) match {

From e0a96e44d0e2a7f7bd97810cc3b6458cb822df9e Mon Sep 17 00:00:00 2001
From: Nathan Fulton <nathan@ibm.com>
Date: Wed, 30 Jun 2021 09:11:45 -0400
Subject: [PATCH 10/14] Adds conversion support for Mathematica's C[].

Mathematica uses C[i] as its default form for parameters/constants.
For example, when computing series approxiamtions to solutions of ODEs, C[i]
is used to represent the constant of integration.

Perhaps a better way of handling this is to ensure that we never get C[i] back from
Mathematica. This can be achieved by always specifying initial conditions in calls to solvers.
For example, instead of:

```
AsymptoticDSolveValue[
  {
    x'[t] == y[t],
    y'[t] == x[t],
    y[0] == x0,
    x[0] == y0
  },
  {x[t], y[t]},
  {t, 0, 10}
]
```

instead do:

```
AsymptoticDSolveValue[
  {
    x'[t] == y[t],
    y'[t] == x[t],
    y[0] == x0,
    x[0] == y0
  },
  {x[t], y[t]},
  {t, 0, 10}
]
```

and then handle x0 and y0 before passing formulas back out of the tool.

Note: this is not currently used in any soundness-critical way that I'm aware of,
but both Reduce and DSolve are capable of generating C[i] symbols. Therefore,
we should do a thorough check that C[i]s are never handed back to us when not expected
or else revert this commit.

See https://reference.wolfram.com/language/ref/C.html for more information on what C[i]
is used for in Mathematica.
---
 .../cmu/cs/ls/keymaerax/tools/qe/MathematicaOpSpec.scala    | 2 ++
 .../cs/ls/keymaerax/tools/qe/MathematicaToKeYmaera.scala    | 6 ++++++
 2 files changed, 8 insertions(+)

diff --git a/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/qe/MathematicaOpSpec.scala b/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/qe/MathematicaOpSpec.scala
index 6751350c94..a7d4719dd5 100644
--- a/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/qe/MathematicaOpSpec.scala
+++ b/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/qe/MathematicaOpSpec.scala
@@ -129,6 +129,8 @@ object MathematicaOpSpec {
 
   def power: BinaryMathOpSpec = BinaryMathOpSpec(symbol("Power"))
 
+  def C: UnaryMathOpSpec = UnaryMathOpSpec(symbol("C")) //@todo document the meaning of this symbol.
+
   // implicit function application name[args]
   def func: NameMathOpSpec = NameMathOpSpec(
     (name: NamedSymbol, args: Array[Expr]) => {
diff --git a/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/qe/MathematicaToKeYmaera.scala b/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/qe/MathematicaToKeYmaera.scala
index cc2a5eb266..eec30bb6b0 100644
--- a/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/qe/MathematicaToKeYmaera.scala
+++ b/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/qe/MathematicaToKeYmaera.scala
@@ -41,6 +41,12 @@ class MathematicaToKeYmaera extends M2KConverter[KExpr] {
     //@note self-created MExpr with head RATIONAL are not rationalQ (type identifiers do not match)
     else if (MathematicaOpSpec.rational.applies(e)) convertBinary(e, Divide.apply)
 
+    // Constant symbols, typically as constants of integration.
+    // Should NOT be leaked out of tooling code. Marking as "interpreted" because that should prevent its unsound use in
+    // at least the most important soundness-critical contexts.
+    // @todo either enforce this invariant or give the function an obviously obnoxious name as a warning to the user.
+    else if (MathematicaOpSpec.C.applies(e))     convertUnary (e, (t: Term) => FuncOf(Function("C", None, Real, Real, true), t))
+
     // Arith expressions
     else if (MathematicaOpSpec.plus.applies(e))   convertNary  (e, Plus.apply)
     else if (MathematicaOpSpec.minus.applies(e))  convertBinary(e, Minus.apply)

From 167c8e120e923220b3d9546ef607d5ba50436a77 Mon Sep 17 00:00:00 2001
From: Nathan Fulton <nathan@ibm.com>
Date: Wed, 30 Jun 2021 09:18:12 -0400
Subject: [PATCH 11/14] Convert DifferentialPrograms to Mathematica.

Note: this code would already written in the DSolve tool, and now
exists in two different places. We should unify this code and clean up
the interface to the DifferenitalProgram converter.
---
 .../qe/DifferentialProgramToMathematica.scala | 81 +++++++++++++++++++
 1 file changed, 81 insertions(+)
 create mode 100644 keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/qe/DifferentialProgramToMathematica.scala

diff --git a/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/qe/DifferentialProgramToMathematica.scala b/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/qe/DifferentialProgramToMathematica.scala
new file mode 100644
index 0000000000..89a0441974
--- /dev/null
+++ b/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/qe/DifferentialProgramToMathematica.scala
@@ -0,0 +1,81 @@
+package edu.cmu.cs.ls.keymaerax.tools.qe
+
+import edu.cmu.cs.ls.keymaerax.core.{AtomicODE, BaseVariable, DifferentialProduct, DifferentialProgram, DifferentialSymbol, Equal, FuncOf, Function, Number, ODESystem, Program, Term, Variable}
+import edu.cmu.cs.ls.keymaerax.infrastruct.ExpressionTraversal.{ExpressionTraversalFunction, StopTraversal}
+import edu.cmu.cs.ls.keymaerax.infrastruct.{ExpressionTraversal, PosInExpr}
+import edu.cmu.cs.ls.keymaerax.tools.ConversionException
+import edu.cmu.cs.ls.keymaerax.tools.ext.ExtMathematicaOpSpec
+import edu.cmu.cs.ls.keymaerax.tools.qe.MathematicaConversion.MExpr
+
+import scala.collection.immutable.{List, Map}
+import scala.math.BigDecimal
+
+case class DifferentialProgramToMathematica(k2m: K2MConverter[MathematicaConversion.KExpr]) {
+  /**
+    * Converts a DifferentialProgram to a Mathematica expression.
+    * This code is largely copy/pasted from MathematicaTools.MathematicalODESolvertool, so that we can KeYmaera
+    * translations of ODEs in other contexts.
+    * @return a 4-tuple containing:
+    *         a list of MExprs specifying the ODEs,
+    *         a list of MExprs specifying symbolic initial conditions,
+    *         a list of MExprs specifying functions (IDK what this is actually -- check the dsolve code and Wolfram documentaiton to figure out)
+    *         the variable used as the time variable (diffArg, as passed in initially)
+    * @author Nathan Fulton
+    */
+  def apply(dp: DifferentialProgram, diffArg: Variable, iv: Map[Variable, Variable]): (List[MExpr], List[MExpr], List[MExpr], Variable) = {
+    /**
+      * @note copied from MathematicaTools.MathematicalODESolvertool.
+      */
+    def toDiffSys(diffSys: DifferentialProgram, diffArg: Variable): List[(Variable, Term)] = {
+      var result = List[(Variable, Term)]()
+      ExpressionTraversal.traverse(new ExpressionTraversalFunction {
+        override def preP(p: PosInExpr, e: Program): Either[Option[StopTraversal], Program] = e match {
+          case AtomicODE(DifferentialSymbol(x), theta) if x != diffArg => result = result :+ (x, theta); Left(None)
+          case AtomicODE(DifferentialSymbol(x), _) if x == diffArg => Left(None)
+          case ODESystem(_, _) => Left(None)
+          case DifferentialProduct(_, _) => Left(None)
+        }
+      }, diffSys)
+      result
+    }
+
+    /** @note coped from dsolve tool. */
+    def functionalizeVars(t: Term, arg: Term, vars: Variable*) = ExpressionTraversal.traverse(
+      new ExpressionTraversalFunction {
+        override def postT(p: PosInExpr, e: Term): Either[Option[StopTraversal], Term] = e match {
+          case v@BaseVariable(name, idx, sort) if vars.isEmpty || vars.contains(v) =>
+            Right(FuncOf(Function(name, idx, arg.sort, sort), arg))
+          case _ => Left(None)
+        }
+      }, t) match {
+      case Some(resultTerm) => resultTerm
+      case None => throw ConversionException("Unable to functionalize " + t)
+    }
+
+    val diffSys = toDiffSys(dp, diffArg)
+
+    val primedVars = diffSys.map(_._1)
+    val functionalizedTerms = diffSys.map{ case (x, theta) => ( x, functionalizeVars(theta, diffArg, primedVars:_*)) }
+    val mathTerms = functionalizedTerms.map({case (x, theta) =>
+      (ExtMathematicaOpSpec.dx(ExtMathematicaOpSpec.primed(k2m(x)))(k2m(diffArg)), k2m(theta))})
+    val convertedDiffSys = mathTerms.map({case (x, theta) => MathematicaOpSpec.equal(x, theta)})
+
+    val functions = diffSys.map(t => k2m(functionalizeVars(t._1, diffArg)))
+
+    //@todo allows for partial initial conditions, but this will probably cause issues if IVs are used.
+    val initialValues = diffSys
+      .map(t =>
+        iv.get(t._1) match {
+          case Some(definedInitialValue) =>
+            Some(
+              Equal(functionalizeVars(t._1, Number(BigDecimal(0)), primedVars:_*), definedInitialValue)
+            )
+          case None => None //@todo perhaps add an error here?
+        }
+      )
+      .filterNot(x => x.isEmpty)
+      .map(x => k2m(x.get))
+
+    (convertedDiffSys, initialValues, functions, diffArg)
+  }
+}

From c20918ca99365143d5092e1df3040988b0fa3341 Mon Sep 17 00:00:00 2001
From: Nathan Fulton <nathan@ibm.com>
Date: Wed, 30 Jun 2021 09:26:18 -0400
Subject: [PATCH 12/14] Tool support for Mathematica's AsymptoticDSolveValue

Computes asymptotic approximations to the solutions of ODEs.

See https://reference.wolfram.com/language/ref/AsymptoticDSolveValue.html
---
 ...ntialSolutionSeriesApproximationTool.scala | 14 ++++
 .../tools/ext/ExtMathematicaOpSpec.scala      |  2 +
 .../ls/keymaerax/tools/ext/Mathematica.scala  | 12 ++-
 .../tools/ext/MathematicaTools.scala          | 75 ++++++++++++++++++-
 .../ls/keymaerax/launcher/TaylorizeMain.scala |  4 +-
 .../tools/SeriesExpansionTests.scala          | 31 ++++++++
 6 files changed, 133 insertions(+), 5 deletions(-)
 create mode 100644 keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/ext/DifferentialSolutionSeriesApproximationTool.scala
 create mode 100644 keymaerax-webui/src/test/scala/edu/cmu/cs/ls/keymaerax/tools/SeriesExpansionTests.scala

diff --git a/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/ext/DifferentialSolutionSeriesApproximationTool.scala b/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/ext/DifferentialSolutionSeriesApproximationTool.scala
new file mode 100644
index 0000000000..0d6697324f
--- /dev/null
+++ b/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/ext/DifferentialSolutionSeriesApproximationTool.scala
@@ -0,0 +1,14 @@
+package edu.cmu.cs.ls.keymaerax.tools.ext
+
+import edu.cmu.cs.ls.keymaerax.core.{NamedSymbol, Number, ODESystem, Term, Variable}
+import edu.cmu.cs.ls.keymaerax.tools.ToolInterface
+
+trait DifferentialSolutionSeriesApproximationTool extends ToolInterface {
+  /**
+    *
+    * @param odes The ODEs
+    * @param ctx Context for any parameters in the ODEs.
+    * @return upper and lower bound series approximations of each primed variable in the ODEs.
+    */
+  def seriesApproximation(odes: ODESystem, ctx: Map[Term, Term]): Option[Map[Variable, Term]]
+}
diff --git a/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/ext/ExtMathematicaOpSpec.scala b/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/ext/ExtMathematicaOpSpec.scala
index 30c4529741..7092a2df97 100644
--- a/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/ext/ExtMathematicaOpSpec.scala
+++ b/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/ext/ExtMathematicaOpSpec.scala
@@ -105,6 +105,8 @@ object ExtMathematicaOpSpec {
 
   def d: BinaryMathOpSpec = BinaryMathOpSpec(symbol("D"))
 
+  def dsolveAsymptoticApproximation: NaryMathOpSpec = NaryMathOpSpec(symbol("AsymptoticDSolveValue"))
+
   //</editor-fold>
 
 }
diff --git a/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/ext/Mathematica.scala b/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/ext/Mathematica.scala
index b96078b070..983c0632e4 100644
--- a/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/ext/Mathematica.scala
+++ b/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/ext/Mathematica.scala
@@ -21,6 +21,7 @@ import edu.cmu.cs.ls.keymaerax.tools.ext.SOSsolveTool.Result
 import edu.cmu.cs.ls.keymaerax.tools.qe.MathematicaConversion.{KExpr, MExpr}
 import edu.cmu.cs.ls.keymaerax.tools.qe.MathematicaOpSpec._
 import edu.cmu.cs.ls.keymaerax.tools.qe.MathematicaToKeYmaera
+import edu.cmu.cs.ls.keymaerax.tools.ext.MathematicaDifferentialSolutionSeriesApproximationTool
 
 import scala.annotation.tailrec
 import scala.collection.immutable.{Map, Seq}
@@ -37,7 +38,7 @@ import scala.collection.mutable.ListBuffer
 class Mathematica(private[tools] val link: MathematicaLink, override val name: String) extends Tool
     with QETacticTool with InvGenTool with ODESolverTool with CounterExampleTool
     with SimulationTool with DerivativeTool with EquationSolverTool with SimplificationTool with AlgebraTool
-    with PDESolverTool with SOSsolveTool with ToolOperationManagement {
+    with PDESolverTool with SOSsolveTool with ToolOperationManagement with DifferentialSolutionSeriesApproximationTool {
 
   /** Indicates whether the tool is initialized. */
   private var initialized = false
@@ -52,6 +53,7 @@ class Mathematica(private[tools] val link: MathematicaLink, override val name: S
   private val mSolve = new MathematicaEquationSolverTool(link)
   private val mAlgebra = new MathematicaAlgebraTool(link)
   private val mSimplify = new MathematicaSimplificationTool(link)
+  private val mExpand = new MathematicaDifferentialSolutionSeriesApproximationTool(link)
 
   private val qeInitialTimeout = Integer.parseInt(Configuration(Configuration.Keys.QE_TIMEOUT_INITIAL))
   private val qeCexTimeout = Integer.parseInt(Configuration(Configuration.Keys.QE_TIMEOUT_CEX))
@@ -70,6 +72,7 @@ class Mathematica(private[tools] val link: MathematicaLink, override val name: S
     mSolve.memoryLimit = memoryLimit
     mAlgebra.memoryLimit = memoryLimit
     mSimplify.memoryLimit = memoryLimit
+    mExpand.memoryLimit = memoryLimit
 
     // initialze tool thread pools
     mQE.init()
@@ -82,6 +85,8 @@ class Mathematica(private[tools] val link: MathematicaLink, override val name: S
     mAlgebra.init()
     mSimplify.init()
     mSOSsolve.init()
+    mExpand.init()
+
 
     initialized = link match {
       case l: JLinkMathematicaLink =>
@@ -113,6 +118,7 @@ class Mathematica(private[tools] val link: MathematicaLink, override val name: S
     mSolve.shutdown()
     mAlgebra.shutdown()
     mSimplify.shutdown()
+    mExpand.shutdown()
     //@note last, because we want to shut down all executors (tool threads) before shutting down the JLink interface
     link match {
       case l: JLinkMathematicaLink => l.shutdown()
@@ -334,4 +340,8 @@ class Mathematica(private[tools] val link: MathematicaLink, override val name: S
 
   /** @inheritdoc */
   override def getAvailableWorkers: Int = 1
+
+  /** @inheritdoc */
+  override def seriesApproximation(odes: ODESystem, ctx: Map[Term, Term]): Option[Map[Variable, Term]] =
+    mExpand.seriesApproximation(odes, ctx)
 }
diff --git a/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/ext/MathematicaTools.scala b/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/ext/MathematicaTools.scala
index 6b9896fe24..127d6f9c13 100644
--- a/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/ext/MathematicaTools.scala
+++ b/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/tools/ext/MathematicaTools.scala
@@ -10,12 +10,13 @@ import edu.cmu.cs.ls.keymaerax.btactics.InvariantGenerator
 import edu.cmu.cs.ls.keymaerax.btactics.helpers.DifferentialHelper
 import edu.cmu.cs.ls.keymaerax.core.{Variable, _}
 import edu.cmu.cs.ls.keymaerax.infrastruct.ExpressionTraversal.{ExpressionTraversalFunction, StopTraversal}
-import edu.cmu.cs.ls.keymaerax.infrastruct.{ExpressionTraversal, FormulaTools, PosInExpr}
+import edu.cmu.cs.ls.keymaerax.infrastruct.{ExpressionTraversal, FormulaTools, PosInExpr, SubstitutionHelper}
 import edu.cmu.cs.ls.keymaerax.tools.qe.MathematicaConversion.{KExpr, _}
 import edu.cmu.cs.ls.keymaerax.tools.ext.SimulationTool.{SimRun, SimState, Simulation}
-import edu.cmu.cs.ls.keymaerax.tools.qe.{BinaryMathOpSpec, ExprFactory, K2MConverter, KeYmaeraToMathematica, M2KConverter, MathematicaNameConversion, MathematicaOpSpec, MathematicaToKeYmaera, NaryMathOpSpec, UnaryMathOpSpec}
+import edu.cmu.cs.ls.keymaerax.tools.qe.{DifferentialProgramToMathematica, BinaryMathOpSpec, ExprFactory, K2MConverter, KeYmaeraToMathematica, M2KConverter, MathematicaNameConversion, MathematicaOpSpec, MathematicaToKeYmaera, NaryMathOpSpec, UnaryMathOpSpec}
 import edu.cmu.cs.ls.keymaerax.tools.qe.MathematicaOpSpec._
 import edu.cmu.cs.ls.keymaerax.tools._
+import edu.cmu.cs.ls.keymaerax.infrastruct.Augmentors._
 
 import scala.collection.immutable
 import scala.math.BigDecimal
@@ -780,6 +781,76 @@ class MathematicaEquationSolverTool(override val link: MathematicaLink) extends
   }
 }
 
+/**
+  * Uses Mathematica's AsymptoticDSolveValue function to construct power series expansions of the solution to a system of ODEs.
+  *
+  * @author Nathan Fulton
+  */
+class MathematicaDifferentialSolutionSeriesApproximationTool(override val link: MathematicaLink) extends BaseKeYmaeraMathematicaBridge[KExpr](link, new UncheckedBaseK2MConverter, PegasusM2KConverter) with DifferentialSolutionSeriesApproximationTool {
+
+  private def pairToList(t: Term): List[Term] = t match {
+    case Pair(left, right) => {
+      assert(!left.isInstanceOf[Pair], "Expected canonical ordering so that tail recursion is possible.")
+      left :: pairToList(right)
+    }
+    case _ => t :: Nil
+  }
+
+  private def convertResult(odes: ODESystem, result: Expression) = {
+    result match {
+      case t: Term => {
+        // Convert nested pairs into a list of terms.
+        val approximations = pairToList(t)
+
+        // Replace all constant symbols in the terms with old(x_i) where x_i is the i^th primed variable in the ODEs.
+        val primedVariables = DifferentialHelper.getPrimedVariables(odes)
+        val approximationsWithInitialConditions = approximations.map(approximation => {
+          Range(1, primedVariables.length + 1).foldLeft(approximation)((approximation, i) => {
+            val constantSymbol = FuncOf(Function("C", None, Real, Real, true), Number(i)) // C[i]
+            val x_i = primedVariables(i - 1) //the i^th primed variable in the system.
+            val old_x_i = FuncOf(Function("old", None, Real, Real), x_i) // old(x_i)
+            approximation.replaceAll(constantSymbol, old_x_i)
+          })
+        })
+
+        Some(
+          DifferentialHelper.getPrimedVariables(odes)
+            .zip(approximationsWithInitialConditions)
+            .toMap
+        )
+      }
+      case _ => throw new ConversionException(s"Expected AsymptoticDSolveValue to give a term or a list of terms, but found ${result}")
+    }
+  }
+
+  /** @inheritdoc */
+  override def seriesApproximation(odes: ODESystem, ctx: Map[Term, Term]): Option[Map[Variable, Term]] = {
+    //See https://reference.wolfram.com/language/ref/AsymptoticDSolveValue.html for the command we're going to use.
+
+    // Apply the context to the ODEs.
+    val odesInCtx = ctx.foldLeft(odes)((currExpr, x) => {
+      SubstitutionHelper.replaceFree(currExpr)(x._1, x._2).asInstanceOf[ODESystem]
+    })
+
+    // Identify a unique time variable.
+    val TIME_VAR = Variable("t") //@todo do something here that's less stupid.
+
+    // Convert the ODEs into Mathematica expressions.
+    // @todo move this out into a utility class independent of the k2m converter interface.
+    val (mODEs, ivs, _, _) = DifferentialProgramToMathematica(k2m).apply(odesInCtx.ode, TIME_VAR, Map())
+
+    //construct the Mathematica expression that will execute AsymptoticDSolveValue.
+    val input = ExtMathematicaOpSpec.dsolveAsymptoticApproximation(
+      MathematicaOpSpec.list(mODEs:_*),
+      MathematicaOpSpec.list(DifferentialHelper.getPrimedVariables(odesInCtx).map(k2m):_*),
+      MathematicaOpSpec.list(k2m(TIME_VAR), k2m(Number(0)), k2m(Number(10))) //@todo I have no clue what these bounds should be. Make these arguments at least.
+    )
+
+    val (_, result) = run(input)
+    convertResult(odes, result)
+  }
+}
+
 /**
   * A link to Mathematica using the JLink interface.
   *
diff --git a/keymaerax-webui/src/main/scala/edu/cmu/cs/ls/keymaerax/launcher/TaylorizeMain.scala b/keymaerax-webui/src/main/scala/edu/cmu/cs/ls/keymaerax/launcher/TaylorizeMain.scala
index 7d29c6a6cd..c542af2436 100644
--- a/keymaerax-webui/src/main/scala/edu/cmu/cs/ls/keymaerax/launcher/TaylorizeMain.scala
+++ b/keymaerax-webui/src/main/scala/edu/cmu/cs/ls/keymaerax/launcher/TaylorizeMain.scala
@@ -11,9 +11,9 @@ import edu.cmu.cs.ls.keymaerax.parser.StringConverter._
 
 
 /**
-  * The main class of Taylorizing sutff from the outside.
+  * The main class for the CLI interface to Taylor approximating solutions to ODEs.
   *
-  * @autor Nathan Fulton
+  * @author Nathan Fulton
   */
 object TaylorizeMain {
   PrettyPrinter.setPrinter(KeYmaeraXPrettyPrinter.pp)
diff --git a/keymaerax-webui/src/test/scala/edu/cmu/cs/ls/keymaerax/tools/SeriesExpansionTests.scala b/keymaerax-webui/src/test/scala/edu/cmu/cs/ls/keymaerax/tools/SeriesExpansionTests.scala
new file mode 100644
index 0000000000..462775ad68
--- /dev/null
+++ b/keymaerax-webui/src/test/scala/edu/cmu/cs/ls/keymaerax/tools/SeriesExpansionTests.scala
@@ -0,0 +1,31 @@
+package edu.cmu.cs.ls.keymaerax.tools
+
+import edu.cmu.cs.ls.keymaerax.btactics.TacticTestBase
+import edu.cmu.cs.ls.keymaerax.core.{NamedSymbol, Number, ODESystem, Term}
+import edu.cmu.cs.ls.keymaerax.parser.StringConverter._
+import smtlib.theories.Core.True
+
+class SeriesExpansionTests extends TacticTestBase {
+  "series expansion tool" should "work for a simple system" in withMathematica(tool => {
+    val odes = ODESystem("{x'=y, y'=-a*x}".asDifferentialProgram, "true".asFormula)
+    val ctx: Map[Term, Term] = Map({"a".asVariable -> Number(5)})
+    val result = tool.seriesApproximation(odes, ctx)
+    println(result)
+  })
+
+  it should "return some result when there's no relevant context." in withMathematica(tool => {
+    val odes = ODESystem("{x'=y, y'=-x}".asDifferentialProgram, "true".asFormula)
+    val ctx: Map[Term, Term] = Map()
+    val result = tool.seriesApproximation(odes, ctx)
+    println(result)
+  })
+
+  it should "return the correct result for x'=1" in withMathematica(tool => {
+    val odes = ODESystem("x'=1".asDifferentialProgram, "true".asFormula)
+    val ctx: Map[Term, Term] = Map()
+    val result = tool.seriesApproximation(odes, ctx)
+
+    result.get("x".asVariable) shouldBe "t + old(x)".asTerm
+  })
+
+}

From f8061968acf374349e4e3e8bc4d8c96fc2daf5ef Mon Sep 17 00:00:00 2001
From: Nathan Fulton <nathan@ibm.com>
Date: Wed, 30 Jun 2021 10:00:26 -0400
Subject: [PATCH 13/14] adds a default tool provider for series approximations.

---
 .../edu/cmu/cs/ls/keymaerax/btactics/ToolProvider.scala   | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/btactics/ToolProvider.scala b/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/btactics/ToolProvider.scala
index 4e86c2255d..0999007545 100644
--- a/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/btactics/ToolProvider.scala
+++ b/keymaerax-core/src/main/scala/edu/cmu/cs/ls/keymaerax/btactics/ToolProvider.scala
@@ -68,6 +68,8 @@ object ToolProvider extends ToolProvider with Logging {
 
   def sosSolveTool(): Option[SOSsolveTool] = f.sosSolveTool()
 
+  def differentialSeriesApproxmationnTool() = f.differentialSeriesApproxmationnTool()
+
   def init(): Boolean = f.init()
 
   def shutdown(): Unit = f.shutdown()
@@ -147,6 +149,9 @@ trait ToolProvider {
   /** Returns a SOSsolve tool. */
   def sosSolveTool(): Option[SOSsolveTool]
 
+  /** Returns a series expansion tool. */
+  def differentialSeriesApproxmationnTool(): Option[DifferentialSolutionSeriesApproximationTool]
+
   /** Initializes the tools. */
   def init(): Boolean
 
@@ -174,6 +179,7 @@ class PreferredToolProvider[T <: Tool](val toolPreferences: List[T]) extends Too
   private[this] lazy val solver: Option[Tool with EquationSolverTool] = toolPreferences.find(_.isInstanceOf[EquationSolverTool]).map(_.asInstanceOf[Tool with EquationSolverTool])
   private[this] lazy val algebra: Option[Tool with AlgebraTool] = toolPreferences.find(_.isInstanceOf[AlgebraTool]).map(_.asInstanceOf[Tool with AlgebraTool])
   private[this] lazy val sossolve: Option[Tool with SOSsolveTool] = toolPreferences.find(_.isInstanceOf[SOSsolveTool]).map(_.asInstanceOf[Tool with SOSsolveTool])
+  private[this] lazy val diffSeriesApproximation: Option[Tool with DifferentialSolutionSeriesApproximationTool] = toolPreferences.find(_.isInstanceOf[DifferentialSolutionSeriesApproximationTool]).map(_.asInstanceOf[Tool with DifferentialSolutionSeriesApproximationTool])
 
   override def tools(): List[Tool] = toolPreferences
   override def defaultTool(): Option[Tool] = toolPreferences.headOption
@@ -193,6 +199,7 @@ class PreferredToolProvider[T <: Tool](val toolPreferences: List[T]) extends Too
   override def solverTool(): Option[EquationSolverTool] = ensureInitialized(solver)
   override def algebraTool(): Option[AlgebraTool] = ensureInitialized(algebra)
   override def sosSolveTool(): Option[SOSsolveTool] = ensureInitialized(sossolve)
+  override def differentialSeriesApproxmationnTool(): Option[DifferentialSolutionSeriesApproximationTool] = ensureInitialized(diffSeriesApproximation)
   override def init(): Boolean = false /* override to initialize tools in more specialized providers */
   override def shutdown(): Unit = toolPreferences.foreach(_.shutdown())
   override def isInitialized: Boolean = toolPreferences.forall(_.isInitialized)
@@ -221,6 +228,7 @@ class NoneToolProvider extends ToolProvider {
   override def solverTool(): Option[EquationSolverTool] = None
   override def algebraTool(): Option[AlgebraTool] = None
   override def sosSolveTool(): Option[SOSsolveTool] = None
+  override def differentialSeriesApproxmationnTool() = None
   override def init(): Boolean = true
   override def shutdown(): Unit = {}
   override def isInitialized: Boolean = true

From 99f066e61f35b7a10d81602a6d4bcb5f0d77a272 Mon Sep 17 00:00:00 2001
From: Nathan Fulton <nathan@ibm.com>
Date: Wed, 30 Jun 2021 10:01:01 -0400
Subject: [PATCH 14/14] incoroproate new series approximation tool into the
 (perhaps poorly named) taylorize cli mode.

---
 .../cmu/cs/ls/keymaerax/launcher/KeYmaeraX.scala |  8 ++++++++
 .../cs/ls/keymaerax/launcher/TaylorizeMain.scala | 16 ++++++++++------
 2 files changed, 18 insertions(+), 6 deletions(-)

diff --git a/keymaerax-webui/src/main/scala/edu/cmu/cs/ls/keymaerax/launcher/KeYmaeraX.scala b/keymaerax-webui/src/main/scala/edu/cmu/cs/ls/keymaerax/launcher/KeYmaeraX.scala
index 6d6a2f6fc6..a59fd9abd0 100644
--- a/keymaerax-webui/src/main/scala/edu/cmu/cs/ls/keymaerax/launcher/KeYmaeraX.scala
+++ b/keymaerax-webui/src/main/scala/edu/cmu/cs/ls/keymaerax/launcher/KeYmaeraX.scala
@@ -97,6 +97,14 @@ object KeYmaeraX {
       //@todo allow multiple passes by filter architecture: -prove bla.key -tactic bla.scal -modelplex -codegen
       options.get('mode) match {
         case Some(Modes.TAYLORIZE) => {
+          val toolConfig =
+            if (options.contains('quantitative)) {
+              configFromFile(Tools.MATHEMATICA) //@note quantitative ModelPlex uses Mathematica to simplify formulas
+            } else {
+              configFromFile("z3")
+            }
+          initializeProver(combineConfigs(options, toolConfig), usage)
+
           val filename = options.get('file)
           filename match {
             case Some(s) => println(TaylorizeMain(s.asInstanceOf[String]))
diff --git a/keymaerax-webui/src/main/scala/edu/cmu/cs/ls/keymaerax/launcher/TaylorizeMain.scala b/keymaerax-webui/src/main/scala/edu/cmu/cs/ls/keymaerax/launcher/TaylorizeMain.scala
index c542af2436..beef72d66f 100644
--- a/keymaerax-webui/src/main/scala/edu/cmu/cs/ls/keymaerax/launcher/TaylorizeMain.scala
+++ b/keymaerax-webui/src/main/scala/edu/cmu/cs/ls/keymaerax/launcher/TaylorizeMain.scala
@@ -4,6 +4,8 @@
 
 package edu.cmu.cs.ls.keymaerax.launcher
 
+import edu.cmu.cs.ls.keymaerax.btactics.ToolProvider
+import edu.cmu.cs.ls.keymaerax.{Configuration, FileConfiguration}
 import edu.cmu.cs.ls.keymaerax.btactics.helpers.DifferentialHelper
 import edu.cmu.cs.ls.keymaerax.core.{DifferentialProgram, ODESystem, PrettyPrinter, Term, Variable}
 import edu.cmu.cs.ls.keymaerax.parser.KeYmaeraXPrettyPrinter
@@ -48,13 +50,15 @@ object TaylorizeMain {
   def taylorize(dp: DifferentialProgram): Option[List[(Variable, Term, Term)]] = {
     val vars = DifferentialHelper.atomicOdes(dp).map(atomic => atomic.xp.x)
 
-    //@todo Grab the taylor approximations from KeYmaera X backend FOR EACH PRIMED VARIABLE.
-    //@todo This whole section is complete nonsesnse filler.
+    val result = ToolProvider.differentialSeriesApproxmationnTool().get.seriesApproximation(
+      ODESystem(dp, "true".asFormula),
+      Map()
+    )
 
-    Some(vars.map((v: Variable) => {
-      val name = v.name
-      (v, s"t^2/2".asTerm, s"t^3/3".asTerm)
-    }))
+    result match {
+      case Some(mapping) => Some(mapping.map(vt => (vt._1,vt._2,vt._2)).toList) //@todo get upper and lower bounds.
+      case None => None
+    }
   }
 
   def output(v: Variable, lowerBound: Term, upperBound: Term) = {