diff --git a/src/compute-engine/differential-equation-utils.ts b/src/compute-engine/differential-equation-utils.ts
index 52ba056a..84e31052 100644
--- a/src/compute-engine/differential-equation-utils.ts
+++ b/src/compute-engine/differential-equation-utils.ts
@@ -1,6 +1,6 @@
 import type { Expression, IComputeEngine } from './global-types';
 import { isFunction, isSymbol, sym } from './boxed-expression/type-guards';
-import { rk4 } from './numerics/differential-equations';
+import { rk4, rk4System } from './numerics/differential-equations';
 
 export function symbolArg(
   engine: IComputeEngine,
@@ -29,34 +29,65 @@ export function isDerivativeOfDependent(
   dependentName: string,
   independentName: string
 ): boolean {
+  return derivativeOrderOfDependent(expr, dependentName, independentName) === 1;
+}
+
+export function derivativeOrderOfDependent(
+  expr: Expression,
+  dependentName: string,
+  independentName: string
+): number | undefined {
   if (isFunction(expr, 'D')) {
-    return (
-      isDependentFunction(expr.op1, dependentName, independentName) &&
-      isSymbol(expr.op2, independentName)
+    const variables = expr.ops.slice(1);
+    if (
+      variables.length === 0 ||
+      !variables.every((op) => isSymbol(op, independentName))
+    )
+      return undefined;
+    const innerOrder = derivativeOrderOfDependent(
+      expr.op1,
+      dependentName,
+      independentName
     );
+    if (innerOrder !== undefined) return innerOrder + variables.length;
+    if (isDependentFunction(expr.op1, dependentName, independentName))
+      return variables.length;
+    return undefined;
   }
 
   if (isFunction(expr, 'Apply') && isFunction(expr.op1, 'Derivative')) {
-    return (
-      isSymbol(expr.op1.op1, dependentName) &&
-      expr.nops === 2 &&
-      isSymbol(expr.op2, independentName)
-    );
+    if (
+      !isSymbol(expr.op1.op1, dependentName) ||
+      expr.nops !== 2 ||
+      !isSymbol(expr.op2, independentName)
+    )
+      return undefined;
+
+    const order = expr.op1.op2 === undefined ? 1 : expr.op1.op2.N().re;
+    return Number.isInteger(order) && order > 0 ? order : undefined;
   }
 
-  return false;
+  return undefined;
 }
 
-function explicitRhs(
+function explicitDerivativeRhs(
   equation: Expression,
   dependentName: string,
   independentName: string
-): Expression | undefined {
+): { order: number; rhs: Expression } | undefined {
   if (!isFunction(equation, 'Equal')) return undefined;
-  if (isDerivativeOfDependent(equation.op1, dependentName, independentName))
-    return equation.op2;
-  if (isDerivativeOfDependent(equation.op2, dependentName, independentName))
-    return equation.op1;
+  const lhsOrder = derivativeOrderOfDependent(
+    equation.op1,
+    dependentName,
+    independentName
+  );
+  if (lhsOrder !== undefined) return { order: lhsOrder, rhs: equation.op2 };
+  const rhsOrder = derivativeOrderOfDependent(
+    equation.op2,
+    dependentName,
+    independentName
+  );
+  if (rhsOrder !== undefined) return { order: rhsOrder, rhs: equation.op1 };
   return undefined;
 }
 
@@ -77,6 +108,32 @@ function substituteDependentCall(
   );
 }
 
+function substituteDependentState(
+  expr: Expression,
+  dependentName: string,
+  independentName: string,
+  stateNames: readonly string[]
+): Expression {
+  if (isDependentFunction(expr, dependentName, independentName))
+    return expr.engine.symbol(stateNames[0]);
+
+  const order = derivativeOrderOfDependent(
+    expr,
+    dependentName,
+    independentName
+  );
+  if (order !== undefined && order < stateNames.length)
+    return expr.engine.symbol(stateNames[order]);
+
+  if (!isFunction(expr)) return expr;
+  return expr.engine._fn(
+    expr.operator,
+    expr.ops.map((op) =>
+      substituteDependentState(op, dependentName, independentName, stateNames)
+    )
+  );
+}
+
 export function nDSolve(
   equation: Expression,
   dependent: Expression,
@@ -92,12 +149,8 @@ export function nDSolve(
   const independentName = sym(limits.op1);
   if (!independentName) return undefined;
 
-  const [x0, x1, y0] = [
-    limits.op2.N().re,
-    limits.op3.N().re,
-    initialValue.N().re,
-  ];
-  if (![x0, x1, y0].every(Number.isFinite)) return undefined;
+  const [x0, x1] = [limits.op2.N().re, limits.op3.N().re];
+  if (![x0, x1].every(Number.isFinite)) return undefined;
 
   const steps = stepsExpr === undefined ? 100 : stepsExpr.N().re;
   if (
@@ -108,12 +161,63 @@ export function nDSolve(
   )
     return undefined;
 
-  const rhs = explicitRhs(equation.structural, dependentName, independentName);
-  if (!rhs) return undefined;
+  const rhsInfo = explicitDerivativeRhs(
+    equation.structural,
+    dependentName,
+    independentName
+  );
+  if (!rhsInfo) return undefined;
+
+  const initialValues = isFunction(initialValue, 'List')
+    ? initialValue.ops.map((op) => op.N().re)
+    : [initialValue.N().re];
+  if (
+    rhsInfo.order !== initialValues.length ||
+    !initialValues.every(Number.isFinite)
+  )
+    return undefined;
+
+  if (rhsInfo.order > 1) {
+    const stateNames = Array.from(
+      { length: rhsInfo.order },
+      (_, i) => `ndsolve${dependentName}state${i}`
+    );
+    const compiledRhs = substituteDependentState(
+      rhsInfo.rhs,
+      dependentName,
+      independentName,
+      stateNames
+    );
+    const compiled = ce._compile(compiledRhs, { realOnly: true });
+    if (!compiled.success) return undefined;
+    const run = compiled.run as (vars: Record<string, number>) => number;
+
+    const samples = rk4System(
+      (x, y) => {
+        const vars: Record<string, number> = { [independentName]: x };
+        stateNames.forEach((name, i) => {
+          vars[name] = y[i];
+        });
+        const highest = run(vars);
+        if (!Number.isFinite(highest)) return undefined;
+        return [...y.slice(1), highest];
+      },
+      x0,
+      initialValues,
+      x1,
+      { steps, deadline: ce._deadline }
+    );
+    if (!samples) return undefined;
+
+    return ce._fn(
+      'List',
+      samples.map(([x, y]) => ce._fn('List', [ce.number(x), ce.number(y[0])]))
+    );
+  }
 
   const stateName = `ndsolve${dependentName}state`;
   const compiledRhs = substituteDependentCall(
-    rhs,
+    rhsInfo.rhs,
     dependentName,
     independentName,
     stateName
@@ -125,7 +229,7 @@ export function nDSolve(
   const samples = rk4(
     (x, y) => run({ [independentName]: x, [stateName]: y }),
     x0,
-    y0,
+    initialValues[0],
     x1,
     { steps, deadline: ce._deadline }
   );
diff --git a/src/compute-engine/numerics/differential-equations.ts b/src/compute-engine/numerics/differential-equations.ts
index 7b5e0f56..bfd6f102 100644
--- a/src/compute-engine/numerics/differential-equations.ts
+++ b/src/compute-engine/numerics/differential-equations.ts
@@ -6,6 +6,7 @@ export type RK4Options = {
 };
 
 export type ODESample = readonly [x: number, y: number];
+export type ODEVectorSample = readonly [x: number, y: readonly number[]];
 
 /**
  * Fixed-step classical fourth-order Runge-Kutta solver for scalar explicit
@@ -50,3 +51,71 @@ export function rk4(
 
   return samples;
 }
+
+/**
+ * Fixed-step classical fourth-order Runge-Kutta solver for first-order systems:
+ * y' = f(x, y), y(x0) = y0.
+ */
+export function rk4System(
+  f: (x: number, y: readonly number[]) => readonly number[] | undefined,
+  x0: number,
+  y0: readonly number[],
+  x1: number,
+  options: RK4Options
+): ODEVectorSample[] | undefined {
+  const steps = Math.trunc(options.steps);
+  if (
+    !Number.isFinite(x0) ||
+    !Number.isFinite(x1) ||
+    !Number.isInteger(steps) ||
+    steps <= 0 ||
+    y0.length === 0 ||
+    !y0.every(Number.isFinite)
+  )
+    return undefined;
+
+  const addScaled = (
+    y: readonly number[],
+    dy: readonly number[],
+    scale: number
+  ): number[] => y.map((yi, i) => yi + scale * dy[i]);
+
+  const combine = (
+    y: readonly number[],
+    k1: readonly number[],
+    k2: readonly number[],
+    k3: readonly number[],
+    k4: readonly number[],
+    h: number
+  ): number[] =>
+    y.map((yi, i) => yi + (h / 6) * (k1[i] + 2 * k2[i] + 2 * k3[i] + k4[i]));
+
+  const h = (x1 - x0) / steps;
+  const samples: ODEVectorSample[] = [[x0, [...y0]]];
+  let x = x0;
+  let y = [...y0];
+
+  for (let i = 0; i < steps; i++) {
+    if ((i & 0xff) === 0) checkDeadline(options.deadline);
+
+    const k1 = f(x, y);
+    if (!k1 || k1.length !== y.length || !k1.every(Number.isFinite))
+      return undefined;
+    const k2 = f(x + h / 2, addScaled(y, k1, h / 2));
+    if (!k2 || k2.length !== y.length || !k2.every(Number.isFinite))
+      return undefined;
+    const k3 = f(x + h / 2, addScaled(y, k2, h / 2));
+    if (!k3 || k3.length !== y.length || !k3.every(Number.isFinite))
+      return undefined;
+    const k4 = f(x + h, addScaled(y, k3, h));
+    if (!k4 || k4.length !== y.length || !k4.every(Number.isFinite))
+      return undefined;
+
+    y = combine(y, k1, k2, k3, k4, h);
+    x = i === steps - 1 ? x1 : x + h;
+    if (!y.every(Number.isFinite)) return undefined;
+    samples.push([x, [...y]]);
+  }
+
+  return samples;
+}
diff --git a/src/compute-engine/symbolic/differential-equations.ts b/src/compute-engine/symbolic/differential-equations.ts
index f2ace10c..83344503 100644
--- a/src/compute-engine/symbolic/differential-equations.ts
+++ b/src/compute-engine/symbolic/differential-equations.ts
@@ -1,10 +1,17 @@
 import { antiderivative } from './antiderivative';
 import type { Expression } from '../global-types';
+import { isValueDef } from '../boxed-expression/utils';
 import { isFunction, isSymbol, sym } from '../boxed-expression/type-guards';
 import {
+  getPolynomialCoefficients,
+  polynomialDegree,
+} from '../boxed-expression/polynomials';
+import {
+  derivativeOrderOfDependent,
   isDependentFunction,
   isDerivativeOfDependent,
 } from '../differential-equation-utils';
+import { durandKernerRoots } from '../numerics/polynomial-roots';
 
 interface LinearTermCoefficients {
   derivative: Expression;
@@ -12,6 +19,11 @@ interface LinearTermCoefficients {
   rest: Expression;
 }
 
+interface DerivativeTermCoefficients {
+  coefficients: Map<number, Expression>;
+  rest: Expression;
+}
+
 function functionName(expr: Expression): string | undefined {
   if (!isFunction(expr)) return undefined;
   return expr.operator;
@@ -80,7 +92,7 @@ function hasDependentOrDerivative(
   independentName: string
 ): boolean {
   if (isDependentFunction(expr, dependentName, independentName)) return true;
-  if (isDerivativeOfDependent(expr, dependentName, independentName))
+  if (derivativeOrderOfDependent(expr, dependentName, independentName))
     return true;
   if (!isFunction(expr)) return false;
   return expr.ops.some((op) =>
@@ -184,22 +196,877 @@ function collectSymbols(
   return symbols;
 }
 
-function integrationConstant(equation: Expression): Expression {
+function freshSymbolName(prefix: string, usedSymbols: Set<string>): string {
+  for (let i = 0; ; i++) {
+    const name = i === 0 ? prefix : `${prefix}_${i}`;
+    if (!usedSymbols.has(name)) return name;
+  }
+}
+
+function integrationConstants(
+  equation: Expression,
+  count: number
+): Expression[] {
   const ce = equation.engine;
   const usedSymbols = collectSymbols(equation);
+  const result: Expression[] = [];
 
-  for (let i = 0; ; i++) {
-    const name = i === 0 ? 'C' : 'c'.repeat(i);
+  for (let i = 1; result.length < count; i++) {
+    const name = `c_${i}`;
     if (usedSymbols.has(name)) continue;
-    if (ce.context.lexicalScope.bindings.has(name)) continue;
-    return ce.symbol(name);
+    const def = ce.lookupDefinition(name);
+    if (def && !(isValueDef(def) && def.value.inferredType)) continue;
+    usedSymbols.add(name);
+    result.push(ce.symbol(name));
+  }
+
+  return result;
+}
+
+function splitDerivativeTerm(
+  term: Expression,
+  dependentName: string,
+  independentName: string
+): { order: number | null; coefficient: Expression } {
+  const ce = term.engine;
+
+  if (isDependentFunction(term, dependentName, independentName))
+    return { order: 0, coefficient: ce.One };
+
+  const directDerivativeOrder = derivativeOrderOfDependent(
+    term,
+    dependentName,
+    independentName
+  );
+  if (directDerivativeOrder !== undefined)
+    return { order: directDerivativeOrder, coefficient: ce.One };
+
+  if (isFunction(term, 'Negate')) {
+    const result = splitDerivativeTerm(
+      term.op1,
+      dependentName,
+      independentName
+    );
+    return { ...result, coefficient: result.coefficient.neg() };
+  }
+
+  if (isFunction(term, 'Multiply')) {
+    let matchedIndex = -1;
+    let matchedOrder: number | null = null;
+
+    for (let i = 0; i < term.ops.length; i++) {
+      let order: number | undefined;
+      if (isDependentFunction(term.ops[i], dependentName, independentName))
+        order = 0;
+      else
+        order = derivativeOrderOfDependent(
+          term.ops[i],
+          dependentName,
+          independentName
+        );
+
+      if (order !== undefined) {
+        if (matchedIndex >= 0) return { order: null, coefficient: term };
+        matchedIndex = i;
+        matchedOrder = order;
+      }
+    }
+
+    if (matchedIndex >= 0 && matchedOrder !== null) {
+      const coefficientFactors = term.ops.filter((_, i) => i !== matchedIndex);
+      const coefficient =
+        coefficientFactors.length === 0
+          ? ce.One
+          : ce.function('Multiply', coefficientFactors);
+      return { order: matchedOrder, coefficient };
+    }
+  }
+
+  return { order: null, coefficient: term };
+}
+
+function collectDerivativeTerms(
+  residual: Expression,
+  dependentName: string,
+  independentName: string
+): DerivativeTermCoefficients {
+  const ce = residual.engine;
+  const terms = isFunction(residual, 'Add')
+    ? residual.ops
+    : isFunction(residual, 'Subtract')
+      ? [residual.op1, residual.op2.neg()]
+      : [residual];
+  const coefficients = new Map<number, Expression>();
+  let rest = ce.Zero;
+
+  for (const term of terms) {
+    const split = splitDerivativeTerm(term, dependentName, independentName);
+    if (split.order === null) rest = rest.add(split.coefficient);
+    else
+      coefficients.set(
+        split.order,
+        (coefficients.get(split.order) ?? ce.Zero).add(split.coefficient)
+      );
+  }
+
+  return { coefficients, rest };
+}
+
+function negateDerivativeCoefficients(
+  coefficients: DerivativeTermCoefficients
+): DerivativeTermCoefficients {
+  const negated = new Map<number, Expression>();
+  for (const [order, coefficient] of coefficients.coefficients)
+    negated.set(order, coefficient.neg());
+  return { coefficients: negated, rest: coefficients.rest.neg() };
+}
+
+function equationDerivativeCoefficients(
+  equation: Expression,
+  dependentName: string,
+  independentName: string
+): DerivativeTermCoefficients {
+  if (!isFunction(equation, 'Equal'))
+    return collectDerivativeTerms(
+      equation.structural,
+      dependentName,
+      independentName
+    );
+
+  const lhs = collectDerivativeTerms(
+    equation.op1.structural,
+    dependentName,
+    independentName
+  );
+  const rhs = negateDerivativeCoefficients(
+    collectDerivativeTerms(
+      equation.op2.structural,
+      dependentName,
+      independentName
+    )
+  );
+  const coefficients = new Map(lhs.coefficients);
+
+  for (const [order, coefficient] of rhs.coefficients)
+    coefficients.set(
+      order,
+      (coefficients.get(order) ?? equation.engine.Zero).add(coefficient)
+    );
+
+  return { coefficients, rest: lhs.rest.add(rhs.rest) };
+}
+
+function isConstantCoefficient(
+  expr: Expression,
+  dependentName: string,
+  independentName: string
+): boolean {
+  if (expr.has(independentName)) return false;
+  return !hasDependentOrDerivative(expr, dependentName, independentName);
+}
+
+function expTerm(coefficient: Expression, root: Expression, variable: string) {
+  const ce = coefficient.engine;
+  const exponent = root.mul(ce.symbol(variable)).simplify();
+  return coefficient.mul(ce.function('Exp', [exponent])).simplify();
+}
+
+function hasOperator(expr: Expression, operator: string): boolean {
+  if (isFunction(expr, operator)) return true;
+  if (!isFunction(expr)) return false;
+  return expr.ops.some((op) => hasOperator(op, operator));
+}
+
+function dSolveAntiderivative(expr: Expression, variable: string): Expression {
+  const ce = expr.engine;
+  if (ce._integrationProvider) {
+    try {
+      const provided = ce._integrationProvider(expr, variable);
+      if (provided && !hasOperator(provided, 'Integrate')) return provided;
+    } catch {
+      // Fall through to the built-in antiderivative, matching Integrate.
+    }
+  }
+  return antiderivative(expr, variable);
+}
+
+function productExpression(ce: Expression['engine'], factors: Expression[]) {
+  if (factors.length === 0) return ce.One;
+  if (factors.length === 1) return factors[0];
+  return ce.function('Multiply', factors).simplify();
+}
+
+function sameArgument(
+  lhs: Expression | undefined,
+  rhs: Expression | undefined
+): boolean {
+  return lhs !== undefined && rhs !== undefined && lhs.isSame(rhs);
+}
+
+function normalizeVariationIntegrand(
+  expr: Expression,
+  independentName: string
+): Expression {
+  const ce = expr.engine;
+  if (isFunction(expr, 'Negate'))
+    return normalizeVariationIntegrand(expr.op1, independentName)
+      .neg()
+      .simplify();
+
+  if (!isFunction(expr, 'Multiply')) return expr;
+
+  const factors = expr.ops.map((op) =>
+    normalizeVariationIntegrand(op, independentName)
+  );
+
+  for (let i = 0; i < factors.length; i++) {
+    for (let j = i + 1; j < factors.length; j++) {
+      const a = factors[i];
+      const b = factors[j];
+      const rest = factors.filter((_, k) => k !== i && k !== j);
+      const trigProduct = normalizedTrigProduct(a, b, independentName);
+      if (!trigProduct) continue;
+      return productExpression(ce, [...rest, trigProduct]).simplify();
+    }
+  }
+
+  return ce.function('Multiply', factors).simplify();
+}
+
+function normalizedTrigProduct(
+  lhs: Expression,
+  rhs: Expression,
+  independentName: string
+): Expression | undefined {
+  const ce = lhs.engine;
+  const ordered =
+    isFunction(lhs, 'Tan') && !isFunction(rhs, 'Tan') ? [rhs, lhs] : [lhs, rhs];
+  const [a, b] = ordered;
+  if (!isFunction(b, 'Tan') || b.op1.has(independentName) === false)
+    return undefined;
+
+  if (isFunction(a, 'Cos') && sameArgument(a.op1, b.op1))
+    return ce.function('Sin', [a.op1]).simplify();
+
+  if (isFunction(a, 'Sin') && sameArgument(a.op1, b.op1)) {
+    const sec = ce.function('Sec', [a.op1]);
+    const cos = ce.function('Cos', [a.op1]);
+    return sec.sub(cos).simplify();
+  }
+
+  return undefined;
+}
+
+function characteristicPolynomial(
+  coefficients: Map<number, Expression>,
+  order: number,
+  variable: string,
+  ce: Expression['engine']
+): Expression {
+  const root = ce.symbol(variable);
+  const terms: Expression[] = [];
+
+  for (let i = 0; i <= order; i++) {
+    const coefficient = (coefficients.get(i) ?? ce.Zero).simplify();
+    if (coefficient.isSame(0)) continue;
+    if (i === 0) terms.push(coefficient);
+    else if (i === 1) terms.push(coefficient.mul(root).simplify());
+    else terms.push(coefficient.mul(root.pow(i)).simplify());
+  }
+
+  return terms.length === 0 ? ce.Zero : ce.function('Add', terms).simplify();
+}
+
+function isZeroAtRoot(
+  expr: Expression,
+  variable: string,
+  root: Expression
+): boolean {
+  const value = expr.subs({ [variable]: root }).simplify();
+  if (value.isSame(0)) return true;
+  const numeric = value.N();
+  return Math.hypot(numeric.re, numeric.im) < 1e-8;
+}
+
+function rootMultiplicity(
+  polynomial: Expression,
+  variable: string,
+  root: Expression,
+  maxMultiplicity: number
+): number {
+  const ce = polynomial.engine;
+  let multiplicity = 0;
+  let derivative = polynomial;
+
+  while (
+    multiplicity < maxMultiplicity &&
+    isZeroAtRoot(derivative, variable, root)
+  ) {
+    multiplicity += 1;
+    derivative = ce
+      .function('D', [derivative, ce.symbol(variable)])
+      .evaluate()
+      .simplify();
+  }
+
+  return multiplicity;
+}
+
+function appendDistinctRoot(
+  roots: Expression[],
+  root: Expression,
+  tolerance = 1e-8
+): void {
+  const rootValue = root.N();
+  if (
+    roots.some((other) => {
+      if (root.isSame(other)) return true;
+      const otherValue = other.N();
+      return (
+        Math.hypot(rootValue.re - otherValue.re, rootValue.im - otherValue.im) <
+        tolerance
+      );
+    })
+  )
+    return;
+
+  roots.push(root);
+}
+
+function numericCharacteristicCoefficients(
+  coefficients: Map<number, Expression>,
+  order: number,
+  ce: Expression['engine']
+): number[] | undefined {
+  const result: number[] = [];
+  for (let i = 0; i <= order; i++) {
+    const value = (coefficients.get(i) ?? ce.Zero).N();
+    if (!Number.isFinite(value.re) || Math.abs(value.im) > 1e-12)
+      return undefined;
+    result.push(value.re);
+  }
+  return result;
+}
+
+function solveHigherOrderWithNumericRoots(
+  equation: Expression,
+  dependentCall: Expression,
+  coefficients: Map<number, Expression>,
+  order: number,
+  independentName: string
+): Expression | undefined {
+  const ce = equation.engine;
+  const coeffs = numericCharacteristicCoefficients(coefficients, order, ce);
+  if (!coeffs) return undefined;
+
+  // CE's exact polynomial root helper can miss complex roots in degree >= 3.
+  // When that happens, use numeric roots rather than returning no solution;
+  // exact recovery of values like sqrt(3)/2 is a future improvement.
+  const roots = durandKernerRoots(coeffs, ce._deadline);
+  if (!roots || roots.length !== order) return undefined;
+
+  const sortedRoots = [...roots].sort((a, b) => {
+    const aReal = Math.abs(a.im) < 1e-8;
+    const bReal = Math.abs(b.im) < 1e-8;
+    if (aReal !== bReal) return aReal ? -1 : 1;
+    if (Math.abs(a.re - b.re) > 1e-8) return a.re - b.re;
+    return a.im - b.im;
+  });
+  const used = Array(sortedRoots.length).fill(false);
+  const constants = integrationConstants(equation, order);
+  const x = ce.symbol(independentName);
+  const terms: Expression[] = [];
+  let constantIndex = 0;
+
+  for (let i = 0; i < sortedRoots.length; i++) {
+    if (used[i]) continue;
+    const root = sortedRoots[i];
+    used[i] = true;
+
+    if (Math.abs(root.im) < 1e-8) {
+      terms.push(
+        expTerm(
+          constants[constantIndex],
+          ce.number(ce.chop(root.re)),
+          independentName
+        )
+      );
+      constantIndex += 1;
+      continue;
+    }
+
+    const conjugateIndex = sortedRoots.findIndex(
+      (candidate, j) =>
+        !used[j] &&
+        Math.abs(candidate.re - root.re) < 1e-7 &&
+        Math.abs(candidate.im + root.im) < 1e-7
+    );
+    if (conjugateIndex < 0 || constantIndex + 1 >= constants.length)
+      return undefined;
+    used[conjugateIndex] = true;
+
+    const alpha = ce.number(ce.chop(root.re));
+    const beta = ce.number(ce.chop(Math.abs(root.im)));
+    const betaX = beta.mul(x).simplify();
+    const oscillatory = constants[constantIndex]
+      .mul(ce.function('Cos', [betaX]))
+      .add(constants[constantIndex + 1].mul(ce.function('Sin', [betaX])))
+      .simplify();
+    terms.push(expTerm(oscillatory, alpha, independentName));
+    constantIndex += 2;
+  }
+
+  if (constantIndex !== order) return undefined;
+  const solution = ce.function('Add', terms).simplify();
+  return ce.function('List', [ce.function('Equal', [dependentCall, solution])]);
+}
+
+function solveHigherOrderHomogeneousConstantCoefficient(
+  equation: Expression,
+  dependentCall: Expression,
+  dependentName: string,
+  independentName: string
+): Expression | undefined {
+  const ce = equation.engine;
+  const collected = equationDerivativeCoefficients(
+    equation,
+    dependentName,
+    independentName
+  );
+  if (!collected.rest.isSame(0)) return undefined;
+
+  const order = Math.max(...collected.coefficients.keys());
+  if (order < 3) return undefined;
+  const leading = (collected.coefficients.get(order) ?? ce.Zero).simplify();
+  if (leading.isSame(0)) return undefined;
+  if (
+    ![...collected.coefficients.values()].every((coefficient) =>
+      isConstantCoefficient(coefficient, dependentName, independentName)
+    )
+  )
+    return undefined;
+
+  const rootVariable = freshSymbolName('dsolveroot', collectSymbols(equation));
+  const polynomial = characteristicPolynomial(
+    collected.coefficients,
+    order,
+    rootVariable,
+    ce
+  );
+  const foundRoots = polynomial.polynomialRoots(rootVariable);
+  if (!foundRoots || foundRoots.length === 0)
+    return solveHigherOrderWithNumericRoots(
+      equation,
+      dependentCall,
+      collected.coefficients,
+      order,
+      independentName
+    );
+
+  const roots: Expression[] = [];
+  for (const root of foundRoots) appendDistinctRoot(roots, root.simplify());
+
+  const rootMultiplicities = roots.map((root) => ({
+    root,
+    multiplicity: rootMultiplicity(polynomial, rootVariable, root, order),
+  }));
+  if (rootMultiplicities.some(({ multiplicity }) => multiplicity === 0))
+    return solveHigherOrderWithNumericRoots(
+      equation,
+      dependentCall,
+      collected.coefficients,
+      order,
+      independentName
+    );
+
+  const totalMultiplicity = rootMultiplicities.reduce(
+    (sum, { multiplicity }) => sum + multiplicity,
+    0
+  );
+  if (totalMultiplicity !== order)
+    return solveHigherOrderWithNumericRoots(
+      equation,
+      dependentCall,
+      collected.coefficients,
+      order,
+      independentName
+    );
+
+  const constants = integrationConstants(equation, order);
+  const x = ce.symbol(independentName);
+  let constantIndex = 0;
+  const terms: Expression[] = [];
+
+  for (const { root, multiplicity } of rootMultiplicities) {
+    for (let power = 0; power < multiplicity; power++) {
+      const coefficient =
+        power === 0
+          ? constants[constantIndex]
+          : constants[constantIndex].mul(x.pow(power)).simplify();
+      terms.push(expTerm(coefficient, root, independentName));
+      constantIndex += 1;
+    }
   }
+
+  const solution = ce.function('Add', terms).simplify();
+  return ce.function('List', [ce.function('Equal', [dependentCall, solution])]);
+}
+
+function solveSecondOrderHomogeneousConstantCoefficient(
+  equation: Expression,
+  dependentCall: Expression,
+  dependentName: string,
+  independentName: string
+): Expression | undefined {
+  const basis = secondOrderConstantCoefficientBasis(
+    equation,
+    dependentName,
+    independentName
+  );
+  if (!basis) return undefined;
+
+  const constants = integrationConstants(equation, 2);
+  const solution = constants[0]
+    .mul(basis[0])
+    .add(constants[1].mul(basis[1]))
+    .simplify();
+  return ceListSolution(dependentCall, solution);
+}
+
+function ceListSolution(
+  dependentCall: Expression,
+  solution: Expression
+): Expression {
+  return dependentCall.engine.function('List', [
+    dependentCall.engine.function('Equal', [dependentCall, solution]),
+  ]);
+}
+
+function secondOrderConstantCoefficientBasis(
+  equation: Expression,
+  dependentName: string,
+  independentName: string
+): [Expression, Expression] | undefined {
+  const collected = equationDerivativeCoefficients(
+    equation,
+    dependentName,
+    independentName
+  );
+  if (!collected.rest.isSame(0)) return undefined;
+  for (const order of collected.coefficients.keys()) {
+    if (order > 2) return undefined;
+  }
+
+  return secondOrderConstantCoefficientBasisFromCoefficients(
+    equation,
+    collected.coefficients,
+    dependentName,
+    independentName
+  );
+}
+
+function secondOrderConstantCoefficientBasisFromCoefficients(
+  equation: Expression,
+  coefficients: Map<number, Expression>,
+  dependentName: string,
+  independentName: string
+): [Expression, Expression] | undefined {
+  const ce = equation.engine;
+  const a = (coefficients.get(2) ?? ce.Zero).simplify();
+  const b = (coefficients.get(1) ?? ce.Zero).simplify();
+  const c0 = (coefficients.get(0) ?? ce.Zero).simplify();
+  if (a.isSame(0)) return undefined;
+  if (
+    ![a, b, c0].every((coefficient) =>
+      isConstantCoefficient(coefficient, dependentName, independentName)
+    )
+  )
+    return undefined;
+
+  const x = ce.symbol(independentName);
+  const twoA = a.mul(2).simplify();
+  const discriminant = b.pow(2).sub(a.mul(c0).mul(4)).simplify();
+
+  if (discriminant.isSame(0)) {
+    const root = b.neg().div(twoA).simplify();
+    const exponential = expTerm(ce.One, root, independentName);
+    return [exponential, x.mul(exponential).simplify()];
+  }
+
+  if (discriminant.isPositive === true) {
+    const sqrtDiscriminant = ce.function('Sqrt', [discriminant]).simplify();
+    const root1 = ce.function('Divide', [
+      ce.function('Add', [b.neg(), sqrtDiscriminant]),
+      twoA,
+    ]);
+    const root2 = ce.function('Divide', [
+      ce.function('Subtract', [b.neg(), sqrtDiscriminant]),
+      twoA,
+    ]);
+    return [
+      expTerm(ce.One, root1, independentName),
+      expTerm(ce.One, root2, independentName),
+    ];
+  }
+
+  if (discriminant.isNegative === true) {
+    const alpha = b.neg().div(twoA).simplify();
+    const beta = ce.function('Sqrt', [discriminant.neg()]).div(twoA).simplify();
+    const exponential = expTerm(ce.One, alpha, independentName);
+    return [
+      exponential.mul(ce.function('Cos', [beta.mul(x).simplify()])).simplify(),
+      exponential.mul(ce.function('Sin', [beta.mul(x).simplify()])).simplify(),
+    ];
+  }
+
+  return undefined;
+}
+
+function solveSecondOrderNonhomogeneousConstantCoefficient(
+  equation: Expression,
+  dependentCall: Expression,
+  dependentName: string,
+  independentName: string
+): Expression | undefined {
+  const ce = equation.engine;
+  const collected = equationDerivativeCoefficients(
+    equation,
+    dependentName,
+    independentName
+  );
+  const a = (collected.coefficients.get(2) ?? ce.Zero).simplify();
+  if (a.isSame(0) || collected.rest.isSame(0)) return undefined;
+  if ([...collected.coefficients.keys()].some((order) => order > 2))
+    return undefined;
+  if (
+    hasDependentOrDerivative(collected.rest, dependentName, independentName) ||
+    ![...collected.coefficients.values()].every((coefficient) =>
+      isConstantCoefficient(coefficient, dependentName, independentName)
+    )
+  )
+    return undefined;
+
+  const basis = secondOrderConstantCoefficientBasisFromCoefficients(
+    equation,
+    collected.coefficients,
+    dependentName,
+    independentName
+  );
+  if (!basis) return undefined;
+
+  const [y1, y2] = basis;
+  const [c1, c2] = integrationConstants(equation, 2);
+  const polynomialParticular = polynomialParticularSolution(
+    equation,
+    collected,
+    independentName
+  );
+  if (polynomialParticular) {
+    const solution = c1
+      .mul(y1)
+      .add(c2.mul(y2))
+      .add(polynomialParticular)
+      .simplify();
+    return ceListSolution(dependentCall, solution);
+  }
+
+  const y1Prime = ce.function('D', [y1, ce.symbol(independentName)]).evaluate();
+  const y2Prime = ce.function('D', [y2, ce.symbol(independentName)]).evaluate();
+  const wronskian = y1.mul(y2Prime).sub(y1Prime.mul(y2)).simplify();
+  if (wronskian.isSame(0)) return undefined;
+
+  const g = collected.rest.neg().simplify();
+  const denominator = a.mul(wronskian).simplify();
+  const u1Integrand = normalizeVariationIntegrand(
+    y2.neg().mul(g).div(denominator).simplify(),
+    independentName
+  );
+  const u2Integrand = normalizeVariationIntegrand(
+    y1.mul(g).div(denominator).simplify(),
+    independentName
+  );
+  const u1 = dSolveAntiderivative(u1Integrand, independentName);
+  const u2 = dSolveAntiderivative(u2Integrand, independentName);
+
+  if (hasOperator(u1, 'Integrate') || hasOperator(u2, 'Integrate'))
+    return undefined;
+
+  const solution = c1
+    .mul(y1)
+    .add(c2.mul(y2))
+    .add(y1.mul(u1))
+    .add(y2.mul(u2))
+    .simplify();
+  return ceListSolution(dependentCall, solution);
+}
+
+function polynomialParticularSolution(
+  equation: Expression,
+  collected: DerivativeTermCoefficients,
+  independentName: string
+): Expression | undefined {
+  const ce = equation.engine;
+  const rhs = collected.rest.neg().simplify();
+  const rhsDegree = polynomialDegree(rhs, independentName);
+  if (rhsDegree < 0) return undefined;
+
+  let zeroRootMultiplicity = 0;
+  while (
+    zeroRootMultiplicity <= 2 &&
+    (collected.coefficients.get(zeroRootMultiplicity) ?? ce.Zero)
+      .simplify()
+      .isSame(0)
+  )
+    zeroRootMultiplicity += 1;
+  if (zeroRootMultiplicity > 2) return undefined;
+
+  const usedSymbols = collectSymbols(equation);
+  const coefficientNames = Array.from({ length: rhsDegree + 1 }, (_, i) => {
+    const name = freshSymbolName(`dsolvep_${i}`, usedSymbols);
+    usedSymbols.add(name);
+    return name;
+  });
+  const x = ce.symbol(independentName);
+  const terms = coefficientNames.map((name, i) =>
+    ce
+      .symbol(name)
+      .mul(x.pow(i + zeroRootMultiplicity))
+      .simplify()
+  );
+  const ansatz = ce.function('Add', terms).simplify();
+
+  const residual = [...collected.coefficients.entries()]
+    .reduce((sum, [order, coefficient]) => {
+      let derivative = ansatz;
+      for (let i = 0; i < order; i++)
+        derivative = ce.function('D', [derivative, x]).evaluate();
+      return sum.add(coefficient.mul(derivative)).simplify();
+    }, collected.rest)
+    .simplify();
+  const residualCoefficients = getPolynomialCoefficients(
+    residual,
+    independentName
+  );
+  if (!residualCoefficients) return undefined;
+
+  const equations = residualCoefficients.map((coefficient) =>
+    ce.function('Equal', [coefficient.simplify(), ce.Zero])
+  );
+  const result = ce.function('List', equations).solve(coefficientNames);
+  const solution = solutionRecord(result);
+  if (!solution) return undefined;
+
+  const particular = ansatz.subs(solution).simplify();
+  if (coefficientNames.some((name) => particular.has(name))) return undefined;
+  return particular;
+}
+
+function solutionRecord(
+  result:
+    | null
+    | ReadonlyArray<Expression>
+    | Record<string, Expression>
+    | Array<Record<string, Expression>>
+): Record<string, Expression> | undefined {
+  if (!result) return undefined;
+  if (!Array.isArray(result)) return result as Record<string, Expression>;
+  const [first] = result;
+  if (!first || 'operator' in first) return undefined;
+  return first as Record<string, Expression>;
+}
+
+function coefficientWithoutPowerOfX(
+  coefficient: Expression,
+  independentName: string,
+  power: number
+): Expression | undefined {
+  const ce = coefficient.engine;
+  const x = ce.symbol(independentName);
+  const xPower = power === 0 ? ce.One : power === 1 ? x : x.pow(power);
+  const scaled = coefficient.div(xPower).simplify();
+  if (scaled.has(independentName)) return undefined;
+  return scaled;
+}
+
+function solveSecondOrderCauchyEulerHomogeneous(
+  equation: Expression,
+  dependentCall: Expression,
+  dependentName: string,
+  independentName: string
+): Expression | undefined {
+  const ce = equation.engine;
+  const collected = equationDerivativeCoefficients(
+    equation,
+    dependentName,
+    independentName
+  );
+  if (!collected.rest.isSame(0)) return undefined;
+  if ([...collected.coefficients.keys()].some((order) => order > 2))
+    return undefined;
+
+  const a = coefficientWithoutPowerOfX(
+    collected.coefficients.get(2) ?? ce.Zero,
+    independentName,
+    2
+  )?.simplify();
+  const b = coefficientWithoutPowerOfX(
+    collected.coefficients.get(1) ?? ce.Zero,
+    independentName,
+    1
+  )?.simplify();
+  const c0 = coefficientWithoutPowerOfX(
+    collected.coefficients.get(0) ?? ce.Zero,
+    independentName,
+    0
+  )?.simplify();
+  if (!a || !b || !c0 || a.isSame(0)) return undefined;
+
+  const [c1, c2] = integrationConstants(equation, 2);
+  const x = ce.symbol(independentName);
+  const effectiveB = b.sub(a).simplify();
+  const twoA = a.mul(2).simplify();
+  const discriminant = effectiveB.pow(2).sub(a.mul(c0).mul(4)).simplify();
+
+  let solution: Expression | undefined;
+  if (discriminant.isSame(0)) {
+    const root = effectiveB.neg().div(twoA).simplify();
+    solution = c1
+      .add(c2.mul(ce.function('Ln', [x])))
+      .mul(x.pow(root))
+      .simplify();
+  } else if (discriminant.isPositive === true) {
+    const sqrtDiscriminant = ce.function('Sqrt', [discriminant]).simplify();
+    const root1 = ce.function('Divide', [
+      ce.function('Add', [effectiveB.neg(), sqrtDiscriminant]),
+      twoA,
+    ]);
+    const root2 = ce.function('Divide', [
+      ce.function('Subtract', [effectiveB.neg(), sqrtDiscriminant]),
+      twoA,
+    ]);
+    solution = c1
+      .mul(x.pow(root1))
+      .add(c2.mul(x.pow(root2)))
+      .simplify();
+  } else if (discriminant.isNegative === true) {
+    const alpha = effectiveB.neg().div(twoA).simplify();
+    const beta = ce.function('Sqrt', [discriminant.neg()]).div(twoA).simplify();
+    const logX = ce.function('Ln', [x]);
+    const oscillatory = c1
+      .mul(ce.function('Cos', [beta.mul(logX).simplify()]))
+      .add(c2.mul(ce.function('Sin', [beta.mul(logX).simplify()])))
+      .simplify();
+    solution = x.pow(alpha).mul(oscillatory).simplify();
+  }
+
+  return solution ? ceListSolution(dependentCall, solution) : undefined;
 }
 
 /**
- * Solve a small first-order linear ODE subset:
+ * Solve a small linear ODE subset:
  *
  *   y'(x) + p(x)y(x) = q(x)
+ *   ay''(x) + by'(x) + cy(x) = 0
  *
  * The returned expression is a `List` of `Equal` expressions for `y(x)`.
  * Unsupported equations return `undefined`, allowing the `DSolve` operator to
@@ -216,6 +1083,41 @@ export function dSolve(
   const { dependentName, independentName, dependentCall } = names;
   const ce = equation.engine;
 
+  const higherOrder = solveHigherOrderHomogeneousConstantCoefficient(
+    equation,
+    dependentCall,
+    dependentName,
+    independentName
+  );
+  if (higherOrder) return higherOrder;
+
+  const cauchyEuler = solveSecondOrderCauchyEulerHomogeneous(
+    equation,
+    dependentCall,
+    dependentName,
+    independentName
+  );
+  if (cauchyEuler) return cauchyEuler;
+
+  const secondOrderNonhomogeneous =
+    solveSecondOrderNonhomogeneousConstantCoefficient(
+      equation,
+      dependentCall,
+      dependentName,
+      independentName
+    );
+  if (secondOrderNonhomogeneous) return secondOrderNonhomogeneous;
+
+  // Keep order 2 separate from the general constant-coefficient solver so
+  // quadratic radical and complex roots stay exact when possible.
+  const secondOrder = solveSecondOrderHomogeneousConstantCoefficient(
+    equation,
+    dependentCall,
+    dependentName,
+    independentName
+  );
+  if (secondOrder) return secondOrder;
+
   const coefficients = equationCoefficients(
     equation,
     dependentName,
@@ -240,17 +1142,17 @@ export function dSolve(
 
   const p = coefficients.dependent.div(coefficients.derivative).simplify();
   const q = coefficients.rest.neg().div(coefficients.derivative).simplify();
-  const c = integrationConstant(equation);
+  const [c] = integrationConstants(equation, 1);
 
   let solution: Expression;
   if (p.isSame(0)) {
-    const integral = antiderivative(q, independentName);
+    const integral = dSolveAntiderivative(q, independentName);
     solution = c.add(integral).simplify();
   } else {
-    const integralP = antiderivative(p, independentName);
+    const integralP = dSolveAntiderivative(p, independentName);
     const integratingFactor = ce.function('Exp', [integralP]).simplify();
     const weightedRhs = integratingFactor.mul(q).simplify();
-    const integral = antiderivative(weightedRhs, independentName);
+    const integral = dSolveAntiderivative(weightedRhs, independentName);
     solution = c.add(integral).div(integratingFactor).simplify();
   }
 
diff --git a/test/compute-engine/differential-equations.test.ts b/test/compute-engine/differential-equations.test.ts
index 2bfea910..69a16fcc 100644
--- a/test/compute-engine/differential-equations.test.ts
+++ b/test/compute-engine/differential-equations.test.ts
@@ -40,7 +40,7 @@ function verifyFirstOrderSolution(
       { match: ['y', 'x'], replace: yValue },
       { recursive: true }
     ) ?? expectedTemplate;
-  const sample = { C: 2, x: 0.75 };
+  const sample = { c_1: 2, x: 0.75 };
   const value = derivative
     .subs(sample)
     .sub(expected.structural.subs(sample))
@@ -49,11 +49,74 @@ function verifyFirstOrderSolution(
   return Math.abs(value) < 1e-10;
 }
 
+function verifyEquationSolution(
+  equation: unknown,
+  solution: ReturnType<typeof dsolve>,
+  sample: Record<string, number>
+): boolean {
+  const solutionEquation = solution.op1;
+  const yValue = solutionEquation.op2;
+  const maxOrder = maxDerivativeOrder(engine.expr(equation, { form: 'raw' }));
+  const derivatives = [yValue];
+  for (let i = 1; i <= maxOrder; i++)
+    derivatives.push(engine.expr(['D', derivatives[i - 1], 'x']).evaluate());
+  const equationExpr = engine.expr(equation, { form: 'raw' });
+  let substituted = equationExpr;
+  for (let order = maxOrder; order >= 1; order--) {
+    let match: unknown = ['y', 'x'];
+    for (let i = 0; i < order; i++) match = ['D', match, 'x'];
+    substituted =
+      substituted.replace(
+        { match, replace: derivatives[order] },
+        { recursive: true }
+      ) ?? substituted;
+    substituted =
+      substituted.replace(
+        {
+          match: ['D', ['y', 'x'], ...Array(order).fill('x')],
+          replace: derivatives[order],
+        },
+        { recursive: true }
+      ) ?? substituted;
+  }
+  substituted =
+    substituted.replace(
+      { match: ['y', 'x'], replace: yValue },
+      { recursive: true }
+    ) ?? substituted;
+  if (substituted.operator !== 'Equal') return false;
+
+  const lhs = substituted.op1.evaluate().canonical;
+  const rhs = substituted.op2.evaluate().canonical;
+  const residual = lhs.sub(rhs).subs(sample).simplify();
+  const value = residual.N().re;
+  return Math.abs(value) < 1e-10;
+}
+
+function maxDerivativeOrder(expr: ReturnType<typeof engine.expr>): number {
+  if (expr.operator === 'D') {
+    let order = expr.ops.length - 1;
+    let inner = expr.op1;
+    while (inner.operator === 'D') {
+      order += inner.ops.length - 1;
+      inner = inner.op1;
+    }
+    return Math.max(order, maxDerivativeOrder(expr.op1));
+  }
+
+  return (expr.ops ?? []).reduce(
+    (max, op) => Math.max(max, maxDerivativeOrder(op)),
+    0
+  );
+}
+
 describe('DSolve', () => {
   test('solves y prime equals y', () => {
     const solution = dsolve(['Equal', ['D', ['y', 'x'], 'x'], ['y', 'x']]);
 
-    expect(solution.toString()).toMatchInlineSnapshot(`[y(x) === C * e^x]`);
+    expect(solution.toString()).toMatchInlineSnapshot(
+      `[y(x) === "c_1" * e^x]`
+    );
     expect(verifyFirstOrderSolution(solution, ['y', 'x'])).toBe(true);
   });
 
@@ -64,7 +127,9 @@ describe('DSolve', () => {
       ['Multiply', 3, ['y', 'x']],
     ]);
 
-    expect(solution.toString()).toMatchInlineSnapshot(`[y(x) === C / e^(-3x)]`);
+    expect(solution.toString()).toMatchInlineSnapshot(
+      `[y(x) === "c_1" / e^(-3x)]`
+    );
     expect(
       verifyFirstOrderSolution(solution, ['Multiply', 3, ['y', 'x']])
     ).toBe(true);
@@ -78,7 +143,7 @@ describe('DSolve', () => {
     ]);
 
     expect(solution.toString()).toMatchInlineSnapshot(
-      `[y(x) === 1/3 * x^3 + C]`
+      `[y(x) === 1/3 * x^3 + "c_1"]`
     );
     expect(verifyFirstOrderSolution(solution, ['Power', 'x', 2])).toBe(true);
   });
@@ -91,7 +156,7 @@ describe('DSolve', () => {
     ]);
 
     expect(solution.toString()).toMatchInlineSnapshot(
-      `[y(x) === x + C / e^x - 1]`
+      `[y(x) === x + "c_1" / e^x - 1]`
     );
     expect(
       verifyFirstOrderSolution(solution, ['Subtract', 'x', ['y', 'x']])
@@ -105,7 +170,9 @@ describe('DSolve', () => {
       0,
     ]);
 
-    expect(solution.toString()).toMatchInlineSnapshot(`[y(x) === C / e^(x^2)]`);
+    expect(solution.toString()).toMatchInlineSnapshot(
+      `[y(x) === "c_1" / e^(x^2)]`
+    );
     expect(
       verifyFirstOrderSolution(solution, [
         'Negate',
@@ -114,15 +181,18 @@ describe('DSolve', () => {
     ).toBe(true);
   });
 
-  test('uses a fallback integration constant when C is already declared', () => {
-    engine.declare('C', 'real');
+  test('uses a fallback integration constant when c_1 is already declared', () => {
+    engine.pushScope();
     try {
+      engine.declare('c_1', 'real');
       const solution = dsolve(['Equal', ['D', ['y', 'x'], 'x'], ['y', 'x']]);
 
-      expect(solution.toString()).toMatchInlineSnapshot(`[y(x) === c * e^x]`);
+      expect(solution.toString()).toMatchInlineSnapshot(
+        `[y(x) === "c_2" * e^x]`
+      );
       expect(verifyFirstOrderSolution(solution, ['y', 'x'])).toBe(true);
     } finally {
-      engine.forget('C');
+      engine.popScope();
     }
   });
 
@@ -136,13 +206,304 @@ describe('DSolve', () => {
     expect(result.operator).toBe('DSolve');
   });
 
-  test('stays inert for unsupported higher-order equations', () => {
+  test('solves second-order homogeneous equation with distinct real roots', () => {
     const result = dsolve([
       'Equal',
       ['D', ['D', ['y', 'x'], 'x'], 'x'],
       ['y', 'x'],
     ]);
 
+    expect(result.toString()).toMatchInlineSnapshot(
+      `[y(x) === "c_1" * e^x + "c_2" * e^(-x)]`
+    );
+    expect(
+      verifyEquationSolution(
+        ['Equal', ['D', ['D', ['y', 'x'], 'x'], 'x'], ['y', 'x']],
+        result,
+        { c_1: 2, c_2: 3, x: 0.75 }
+      )
+    ).toBe(true);
+  });
+
+  test('solves flat-form second-order derivatives', () => {
+    const equation = ['Equal', ['D', ['y', 'x'], 'x', 'x'], ['y', 'x']];
+    const result = dsolve(equation);
+
+    expect(result.toString()).toMatchInlineSnapshot(
+      `[y(x) === "c_1" * e^x + "c_2" * e^(-x)]`
+    );
+    expect(
+      verifyEquationSolution(equation, result, { c_1: 2, c_2: 3, x: 0.75 })
+    ).toBe(true);
+  });
+
+  test('keeps irrational characteristic roots exact when possible', () => {
+    const equation = [
+      'Equal',
+      [
+        'Add',
+        ['D', ['D', ['y', 'x'], 'x'], 'x'],
+        ['Negate', ['D', ['y', 'x'], 'x']],
+        ['Negate', ['y', 'x']],
+      ],
+      0,
+    ];
+    const solution = dsolve(equation);
+
+    expect(solution.toString()).toMatchInlineSnapshot(
+      `[y(x) === "c_1" * e^(1/2 * x * (1 + sqrt(5))) + "c_2" * e^(1/2 * x * (1 - sqrt(5)))]`
+    );
+    expect(
+      verifyEquationSolution(equation, solution, { c_1: 2, c_2: 3, x: 0.75 })
+    ).toBe(true);
+  });
+
+  test('solves second-order homogeneous equation with complex roots', () => {
+    const equation = [
+      'Equal',
+      ['Add', ['D', ['D', ['y', 'x'], 'x'], 'x'], ['y', 'x']],
+      0,
+    ];
+    const solution = dsolve(equation);
+
+    expect(solution.toString()).toMatchInlineSnapshot(
+      `[y(x) === "c_1" * cos(x) + "c_2" * sin(x)]`
+    );
+    expect(
+      verifyEquationSolution(equation, solution, { c_1: 2, c_2: 3, x: 0.75 })
+    ).toBe(true);
+  });
+
+  test('solves second-order homogeneous equation with a repeated root', () => {
+    const equation = [
+      'Equal',
+      [
+        'Add',
+        ['D', ['D', ['y', 'x'], 'x'], 'x'],
+        ['Negate', ['Multiply', 2, ['D', ['y', 'x'], 'x']]],
+        ['y', 'x'],
+      ],
+      0,
+    ];
+    const solution = dsolve(equation);
+
+    expect(solution.toString()).toMatchInlineSnapshot(
+      `[y(x) === "c_2" * x * e^x + "c_1" * e^x]`
+    );
+    expect(
+      verifyEquationSolution(equation, solution, { c_1: 2, c_2: 3, x: 0.75 })
+    ).toBe(true);
+  });
+
+  test('solves second-order homogeneous equation with zero repeated root', () => {
+    const equation = ['Equal', ['D', ['D', ['y', 'x'], 'x'], 'x'], 0];
+    const solution = dsolve(equation);
+
+    expect(solution.toString()).toMatchInlineSnapshot(
+      `[y(x) === "c_2" * x + "c_1"]`
+    );
+    expect(
+      verifyEquationSolution(equation, solution, { c_1: 2, c_2: 3, x: 0.75 })
+    ).toBe(true);
+  });
+
+  test('solves third-order homogeneous equation with real roots', () => {
+    const equation = [
+      'Equal',
+      [
+        'Add',
+        ['D', ['D', ['D', ['y', 'x'], 'x'], 'x'], 'x'],
+        ['Negate', ['Multiply', 6, ['D', ['D', ['y', 'x'], 'x'], 'x']]],
+        ['Multiply', 11, ['D', ['y', 'x'], 'x']],
+        ['Negate', ['Multiply', 6, ['y', 'x']]],
+      ],
+      0,
+    ];
+    const solution = dsolve(equation);
+
+    expect(solution.toString()).toMatchInlineSnapshot(
+      `[y(x) === "c_1" * e^x + "c_2" * e^(2x) + "c_3" * e^(3x)]`
+    );
+    expect(
+      verifyEquationSolution(equation, solution, {
+        c_1: 2,
+        c_2: 3,
+        c_3: 5,
+        x: 0.25,
+      })
+    ).toBe(true);
+  });
+
+  test('solves third-order homogeneous equation with repeated root', () => {
+    const equation = [
+      'Equal',
+      [
+        'Add',
+        ['D', ['D', ['D', ['y', 'x'], 'x'], 'x'], 'x'],
+        ['Negate', ['Multiply', 3, ['D', ['D', ['y', 'x'], 'x'], 'x']]],
+        ['Multiply', 3, ['D', ['y', 'x'], 'x']],
+        ['Negate', ['y', 'x']],
+      ],
+      0,
+    ];
+    const solution = dsolve(equation);
+
+    expect(solution.toString()).toMatchInlineSnapshot(
+      `[y(x) === "c_3" * x^2 * e^x + "c_2" * x * e^x + "c_1" * e^x]`
+    );
+    expect(
+      verifyEquationSolution(equation, solution, {
+        c_1: 2,
+        c_2: 3,
+        c_3: 5,
+        x: 0.25,
+      })
+    ).toBe(true);
+  });
+
+  test('solves third-order homogeneous equation with numeric complex roots', () => {
+    const equation = [
+      'Equal',
+      [
+        'Add',
+        ['D', ['D', ['D', ['y', 'x'], 'x'], 'x'], 'x'],
+        ['Negate', ['y', 'x']],
+      ],
+      0,
+    ];
+    const solution = dsolve(equation);
+
+    expect(solution.toString()).toMatchInlineSnapshot(
+      `[y(x) === "c_1" * e^x + "c_2" * cos(0.8660254037844387 * x) * e^(-0.5 * x) + "c_3" * sin(0.8660254037844387 * x) * e^(-0.5 * x)]`
+    );
+    expect(
+      verifyEquationSolution(equation, solution, {
+        c_1: 2,
+        c_2: 3,
+        c_3: 5,
+        x: 0.25,
+      })
+    ).toBe(true);
+  });
+
+  test('solves nonhomogeneous second-order constant coefficient equation', () => {
+    const equation = [
+      'Equal',
+      ['D', ['D', ['y', 'x'], 'x'], 'x'],
+      1,
+    ];
+    const result = dsolve(equation);
+
+    expect(result.toString()).toMatchInlineSnapshot(
+      `[y(x) === 1/2 * x^2 + "c_2" * x + "c_1"]`
+    );
+    expect(
+      verifyEquationSolution(equation, result, { c_1: 2, c_2: 3, x: 0.75 })
+    ).toBe(true);
+  });
+
+  test('solves polynomial-forced second-order constant coefficient equation', () => {
+    const equation = [
+      'Equal',
+      [
+        'Add',
+        ['D', ['D', ['y', 'x'], 'x'], 'x'],
+        ['Negate', ['y', 'x']],
+      ],
+      'x',
+    ];
+    const result = dsolve(equation);
+
+    expect(result.toString()).toMatchInlineSnapshot(
+      `[y(x) === -x + "c_1" * e^x + "c_2" * e^(-x)]`
+    );
+    expect(
+      verifyEquationSolution(equation, result, { c_1: 2, c_2: 3, x: 0.75 })
+    ).toBe(true);
+  });
+
+  test('solves tangent-forced second-order constant coefficient equation', () => {
+    const equation = [
+      'Equal',
+      ['Add', ['D', ['D', ['y', 'x'], 'x'], 'x'], ['y', 'x']],
+      ['Tan', 'x'],
+    ];
+    const result = dsolve(equation);
+
+    expect(result.toString()).toMatchInlineSnapshot(
+      `[y(x) === "c_1" * cos(x) + "c_2" * sin(x) - cos(x) * ln(|tan(x) + sec(x)|)]`
+    );
+    expect(result.operator).toBe('List');
+  });
+
+  test('solves Cauchy-Euler equation with distinct roots', () => {
+    const equation = [
+      'Equal',
+      [
+        'Add',
+        ['Multiply', ['Power', 'x', 2], ['D', ['D', ['y', 'x'], 'x'], 'x']],
+        ['Negate', ['Multiply', 2, ['y', 'x']]],
+      ],
+      0,
+    ];
+    const result = dsolve(equation);
+
+    expect(result.toString()).toMatchInlineSnapshot(
+      `[y(x) === "c_1" * x^2 + "c_2" / x]`
+    );
+    expect(
+      verifyEquationSolution(equation, result, { c_1: 2, c_2: 3, x: 2 })
+    ).toBe(true);
+  });
+
+  test('solves Cauchy-Euler equation with repeated roots', () => {
+    const equation = [
+      'Equal',
+      [
+        'Add',
+        ['Multiply', ['Power', 'x', 2], ['D', ['D', ['y', 'x'], 'x'], 'x']],
+        ['Multiply', 'x', ['D', ['y', 'x'], 'x']],
+      ],
+      0,
+    ];
+    const result = dsolve(equation);
+
+    expect(result.toString()).toMatchInlineSnapshot(
+      `[y(x) === "c_1" + "c_2" * ln(x)]`
+    );
+    expect(
+      verifyEquationSolution(equation, result, { c_1: 2, c_2: 3, x: 2 })
+    ).toBe(true);
+  });
+
+  test('solves Cauchy-Euler equation with complex roots', () => {
+    const equation = [
+      'Equal',
+      [
+        'Add',
+        ['Multiply', ['Power', 'x', 2], ['D', ['D', ['y', 'x'], 'x'], 'x']],
+        ['Multiply', 'x', ['D', ['y', 'x'], 'x']],
+        ['y', 'x'],
+      ],
+      0,
+    ];
+    const result = dsolve(equation);
+
+    expect(result.toString()).toMatchInlineSnapshot(
+      `[y(x) === "c_1" * cos(ln(x)) + "c_2" * sin(ln(x))]`
+    );
+    expect(
+      verifyEquationSolution(equation, result, { c_1: 2, c_2: 3, x: 2 })
+    ).toBe(true);
+  });
+
+  test('stays inert when variation of parameters cannot integrate', () => {
+    const result = dsolve([
+      'Equal',
+      ['Add', ['D', ['D', ['y', 'x'], 'x'], 'x'], ['y', 'x']],
+      ['Divide', 1, 'x'],
+    ]);
+
     expect(result.operator).toBe('DSolve');
   });
 });
@@ -217,6 +578,40 @@ describe('NDSolve', () => {
     expect(y).toBeCloseTo(expected, 10);
   });
 
+  test('solves second-order IVP with RK4 system samples', () => {
+    const result = ndsolve(
+      ['Equal', ['D', ['D', ['y', 'x'], 'x'], 'x'], ['Negate', ['y', 'x']]],
+      ['List', 0, 1],
+      200
+    );
+    const [, y] = finalSample(result);
+
+    expect(result.operator).toBe('List');
+    expect(y).toBeCloseTo(Math.sin(1), 10);
+  });
+
+  test('solves third-order IVP with RK4 system samples', () => {
+    const result = ndsolve(
+      ['Equal', ['D', ['D', ['D', ['y', 'x'], 'x'], 'x'], 'x'], ['y', 'x']],
+      ['List', 1, 1, 1],
+      200
+    );
+    const [, y] = finalSample(result);
+
+    expect(result.operator).toBe('List');
+    expect(y).toBeCloseTo(Math.E, 10);
+  });
+
+  test('stays inert when higher-order IVP initial values have wrong length', () => {
+    const result = ndsolve(
+      ['Equal', ['D', ['D', ['y', 'x'], 'x'], 'x'], ['Negate', ['y', 'x']]],
+      ['List', 0],
+      200
+    );
+
+    expect(result.operator).toBe('NDSolve');
+  });
+
   test('stays inert for unsupported implicit equations', () => {
     const result = ndsolve(
       ['Equal', ['Add', ['D', ['y', 'x'], 'x'], ['y', 'x']], 'x'],