diff --git a/src/include/serif/utils/laneemden/laneemden.hpp b/src/include/serif/utils/laneemden/laneemden.hpp
new file mode 100644
index 00000000..b339a097
--- /dev/null
+++ b/src/include/serif/utils/laneemden/laneemden.hpp
@@ -0,0 +1,4 @@
+#pragma once
+
+#include "serif/utils/laneemden/solution.hpp"
+#include "serif/utils/laneemden/scaling.hpp"
diff --git a/src/include/serif/utils/laneemden/scaling.hpp b/src/include/serif/utils/laneemden/scaling.hpp
new file mode 100644
index 00000000..8a25ee52
--- /dev/null
+++ b/src/include/serif/utils/laneemden/scaling.hpp
@@ -0,0 +1,63 @@
+#pragma once
+
+#include <expected>
+#include <optional>
+#include <string>
+
+#include "serif/dimensions/dimensions.hpp"
+#include "serif/utils/laneemden/solution.hpp"
+
+namespace serif::utils::laneemden {
+    /** @brief A physical conversion of one finite-radius dimensionless polytrope, in the caller's consistent units. */
+    struct PhysicalScales final {
+        double n;
+        double radial_scale;
+        dimensions::QuantityValue<dimensions::Mass> mass;
+        dimensions::QuantityValue<dimensions::Length> radius;
+        dimensions::QuantityValue<dimensions::MassDensity> central_density;
+        dimensions::QuantityValue<dimensions::SpecificEnthalpy> central_enthalpy;
+        dimensions::QuantityValue<dimensions::GravitationalConstant> gravitational_constant;
+        // K has n-dependent dimensions; absent for n = 0 or when this optional derived value is not representable.
+        std::optional<double> polytropic_constant;
+    };
+
+    /** @brief Physical scalar fields; gravity is the outward potential derivative, not the inward force component. */
+    struct PhysicalSample final {
+        dimensions::QuantityValue<dimensions::MassDensity> density;
+        dimensions::QuantityValue<dimensions::Pressure> pressure;
+        dimensions::QuantityValue<dimensions::SpecificEnthalpy> enthalpy;
+        dimensions::QuantityValue<dimensions::GravitationalPotential> potential;
+        dimensions::QuantityValue<dimensions::Mass> enclosed_mass;
+        dimensions::QuantityValue<dimensions::Acceleration> gravity;
+    };
+
+    /** @brief Fix a reference star by mass and central density; its implied K is an output, never a change to an EOS. */
+    [[nodiscard]] std::expected<PhysicalScales, std::string> scale_from_mass(
+        const Solution &solution,
+        dimensions::QuantityValue<dimensions::Mass> mass,
+        dimensions::QuantityValue<dimensions::MassDensity> central_density,
+        dimensions::QuantityValue<dimensions::GravitationalConstant> gravitational_constant);
+
+    /** @brief Fix a physical star by central enthalpy and density, including the incompressible n = 0 limit. */
+    [[nodiscard]] std::expected<PhysicalScales, std::string> scale_from_enthalpy(
+        const Solution &solution,
+        dimensions::QuantityValue<dimensions::SpecificEnthalpy> central_enthalpy,
+        dimensions::QuantityValue<dimensions::MassDensity> central_density,
+        dimensions::QuantityValue<dimensions::GravitationalConstant> gravitational_constant);
+
+    /** @brief Fix a polytropic EOS through its K and central density; n must be strictly positive. */
+    [[nodiscard]] std::expected<PhysicalScales, std::string> scale_from_polytropic_constant(
+        const Solution &solution,
+        double polytropic_constant,
+        dimensions::QuantityValue<dimensions::MassDensity> central_density,
+        dimensions::QuantityValue<dimensions::GravitationalConstant> gravitational_constant);
+
+    /**
+     * @brief Evaluate physical fields, continuing into vacuum beyond the stellar surface.
+     *
+     * The potential vanishes at infinity. Density, pressure, and enthalpy are zero at and
+     * beyond the surface, including the exterior side of the incompressible density jump.
+     */
+    [[nodiscard]] std::expected<PhysicalSample, std::string> evaluate_physical(
+        const Solution &solution, const PhysicalScales &scales, dimensions::QuantityValue<dimensions::Length> radius);
+}
diff --git a/src/include/serif/utils/laneemden/solution.hpp b/src/include/serif/utils/laneemden/solution.hpp
new file mode 100644
index 00000000..096aa168
--- /dev/null
+++ b/src/include/serif/utils/laneemden/solution.hpp
@@ -0,0 +1,141 @@
+#pragma once
+
+#include <cstddef>
+#include <expected>
+#include <optional>
+#include <string>
+#include <string_view>
+#include <vector>
+#include <format>
+
+#include "serif/utils/error/errors.hpp"
+
+namespace serif::utils::laneemden {
+    /** @brief Dimensionless integration controls; the surface tolerance does not bound the global integration error. */
+    struct IntegrationOptions final {
+        double step = 0.01;
+        double maximum_radius = 100.0;
+        std::size_t maximum_steps = 100000;
+        double surface_tolerance = 1e-10;
+        std::size_t central_series_terms = 12;
+
+        /** @brief Compare the numerical prescription when tracking reference dependencies. */
+        friend bool operator==(const IntegrationOptions &, const IntegrationOptions &) = default;
+    };
+
+    /** @brief The regular Lane–Emden solution and its first derivative at one dimensionless radius. */
+    struct Sample final {
+        double xi;
+        double theta;
+        double derivative;
+    };
+
+    /** @brief A regular profile, ending at its first zero or the requested integration radius. */
+    struct Solution final {
+        double n;
+        std::vector<Sample> samples;
+        std::optional<Sample> surface;
+    };
+
+    /**
+     * @brief Integrate theta'' + 2 theta'/xi + theta^n = 0 from the regular center.
+     *
+     * The central series follows Boudreaux's original power-series recurrence. Both dependent
+     * variables now advance through the same four Runge–Kutta stages. A missing surface is a
+     * valid finite-interval result, including the n >= 5 solutions. Exhausting the step budget
+     * or producing an unrepresentable numerical state returns a failure instead.
+     */
+    [[nodiscard]] std::expected<Solution, std::string> solve(double n, IntegrationOptions options = {});
+
+    /** @brief Evaluate a solved dimensionless profile with cubic Hermite interpolation, without extrapolation. */
+    [[nodiscard]] std::expected<Sample, std::string> evaluate(const Solution &solution, double xi);
+}
+
+template <>
+struct std::formatter<serif::utils::laneemden::IntegrationOptions> {
+    template <typename FormatContext>
+    auto format(const serif::utils::laneemden::IntegrationOptions &options, FormatContext &ctx) const {
+        return std::format_to(ctx.out(), "step={}, maximum_radius={}, maximum_steps={}, surface_tolerance={}, central_series_terms={}",
+            options.step, options.maximum_radius, options.maximum_steps, options.surface_tolerance, options.central_series_terms);
+    }
+};
+
+/** @brief Select all sample values, or one named component with :xi, :theta, or :derivative. */
+template <>
+struct std::formatter<serif::utils::laneemden::Sample> {
+    enum class Component { all, xi, theta, derivative };
+    Component component{Component::all};
+
+    /** @brief Validate the selector during literal format checking or runtime format parsing. */
+    constexpr auto parse(std::format_parse_context &ctx) {
+        using namespace serif::utils::errors;
+
+        component = Component::all;
+        auto end = ctx.begin();
+        while (end != ctx.end() && *end != '}') {
+            ++end;
+        }
+        if (end == ctx.begin()) {
+            return end;
+        }
+        const std::string_view selection(ctx.begin(), end);
+        if (selection == "all") {
+            component = Component::all;
+        } else if (selection == "xi") {
+            component = Component::xi;
+        } else if (selection == "theta") {
+            component = Component::theta;
+        } else if (selection == "derivative") {
+            component = Component::derivative;
+        } else {
+            serif_error<std::format_error, SERiFErrorCode::INVALID_ARGUMENT>(
+                "Unknown Lane–Emden sample selector '{}'. Use an empty specifier, 'all', 'xi', 'theta', or 'derivative'.", selection);
+        }
+        return end;
+    }
+
+    /** @brief Keep the labeled tuple for complete samples and emit a bare number for a selected component. */
+    template <typename FormatContext>
+    auto format(const serif::utils::laneemden::Sample &sample, FormatContext &ctx) const {
+        switch (component) {
+            case Component::xi:
+                return std::format_to(ctx.out(), "{}", sample.xi);
+            case Component::theta:
+                return std::format_to(ctx.out(), "{}", sample.theta);
+            case Component::derivative:
+                return std::format_to(ctx.out(), "{}", sample.derivative);
+            default:
+                return std::format_to(ctx.out(), "(ξ={}, θ={}, dθ/dξ={})", sample.xi, sample.theta, sample.derivative);
+        }
+    }
+};
+
+/** @brief Describe the solution while applying one shared selection to its profile and optional surface sample. */
+template <>
+struct std::formatter<serif::utils::laneemden::Solution> {
+    std::formatter<serif::utils::laneemden::Sample> sample_formatter;
+
+    /** @brief Delegate selection to Sample so the two formatters cannot develop different grammars. */
+    constexpr auto parse(std::format_parse_context &ctx) {
+        return sample_formatter.parse(ctx);
+    }
+
+    /** @brief Stream the profile through its parsed sample formatter without building intermediate strings. */
+    template <typename FormatContext>
+    auto format(const serif::utils::laneemden::Solution &solution, FormatContext &ctx) const {
+        ctx.advance_to(std::format_to(ctx.out(), "Solution(n={}, samples=[", solution.n));
+        for (std::size_t index = 0; index < solution.samples.size(); ++index) {
+            if (index != 0) {
+                ctx.advance_to(std::format_to(ctx.out(), ", "));
+            }
+            ctx.advance_to(sample_formatter.format(solution.samples[index], ctx));
+        }
+        ctx.advance_to(std::format_to(ctx.out(), "], surface="));
+        if (solution.surface) {
+            ctx.advance_to(sample_formatter.format(*solution.surface, ctx));
+        } else {
+            ctx.advance_to(std::format_to(ctx.out(), "none"));
+        }
+        return std::format_to(ctx.out(), ")");
+    }
+};
diff --git a/src/lib/serif/utils/laneemden/scaling.cpp b/src/lib/serif/utils/laneemden/scaling.cpp
new file mode 100644
index 00000000..9ea762c7
--- /dev/null
+++ b/src/lib/serif/utils/laneemden/scaling.cpp
@@ -0,0 +1,139 @@
+#include "serif/utils/laneemden/scaling.hpp"
+
+#include <algorithm>
+#include <cmath>
+#include <numbers>
+
+namespace serif::utils::laneemden {
+    namespace {
+        constexpr long double four_pi = 4.0L * std::numbers::pi_v<long double>;
+
+        /** @brief Validate a positive physical scale without assuming a particular unit system. */
+        bool positive_finite(double value) {
+            return std::isfinite(value) && value > 0.0;
+        }
+
+        /** @brief Require the finite surface needed to assign a stellar radius and total mass. */
+        bool finite_surface(const Solution &solution) {
+            return std::isfinite(solution.n) && solution.n >= 0.0 && solution.surface && positive_finite(solution.surface->xi) &&
+                std::isfinite(solution.surface->derivative) && solution.surface->derivative < 0.0;
+        }
+
+        /** @brief Collect the shared scaling identities, keeping intermediate products representable in physical units. */
+        std::expected<PhysicalScales, std::string> make_scales(
+            const Solution &solution, long double radial_scale, double central_density, double gravitational_constant) {
+            const long double surface_xi = solution.surface->xi;
+            const long double dimensionless_mass = -surface_xi * surface_xi * solution.surface->derivative;
+            const double radius = static_cast<double>(radial_scale * surface_xi);
+            const double mass = static_cast<double>(four_pi * central_density * radial_scale * radial_scale * radial_scale * dimensionless_mass);
+            const double enthalpy = static_cast<double>(four_pi * gravitational_constant * central_density * radial_scale * radial_scale);
+            const double alpha = static_cast<double>(radial_scale);
+            if (!positive_finite(alpha) || !positive_finite(radius) || !positive_finite(mass) || !positive_finite(enthalpy)) {
+                return std::unexpected("Lane–Emden physical scales overflow or underflow; choose representable consistent units.");
+            }
+            std::optional<double> polytropic_constant;
+            if (solution.n > 0.0) {
+                const double constant = static_cast<double>(enthalpy / ((solution.n + 1.0L) * std::pow(static_cast<long double>(central_density), 1.0L / solution.n)));
+                if (positive_finite(constant)) {
+                    polytropic_constant = constant;
+                }
+            }
+            return PhysicalScales{
+                .n = solution.n,
+                .radial_scale = alpha,
+                .mass = dimensions::QuantityValue<dimensions::Mass>{mass},
+                .radius = dimensions::QuantityValue<dimensions::Length>{radius},
+                .central_density = dimensions::QuantityValue<dimensions::MassDensity>{central_density},
+                .central_enthalpy = dimensions::QuantityValue<dimensions::SpecificEnthalpy>{enthalpy},
+                .gravitational_constant = dimensions::QuantityValue<dimensions::GravitationalConstant>{gravitational_constant},
+                .polytropic_constant = polytropic_constant
+            };
+        }
+    }
+
+    std::expected<PhysicalScales, std::string> scale_from_mass(
+        const Solution &solution,
+        dimensions::QuantityValue<dimensions::Mass> mass,
+        dimensions::QuantityValue<dimensions::MassDensity> central_density,
+        dimensions::QuantityValue<dimensions::GravitationalConstant> gravitational_constant) {
+        if (!finite_surface(solution) || !positive_finite(mass.value()) || !positive_finite(central_density.value()) || !positive_finite(gravitational_constant.value())) {
+            return std::unexpected("Lane–Emden mass scaling needs a finite surface and positive finite mass, central density, and G.");
+        }
+        const long double xi = solution.surface->xi;
+        const long double dimensionless_mass = -xi * xi * solution.surface->derivative;
+        const long double alpha = std::cbrt(mass.value() / (four_pi * central_density.value() * dimensionless_mass));
+        return make_scales(solution, alpha, central_density.value(), gravitational_constant.value());
+    }
+
+    std::expected<PhysicalScales, std::string> scale_from_enthalpy(
+        const Solution &solution,
+        dimensions::QuantityValue<dimensions::SpecificEnthalpy> central_enthalpy,
+        dimensions::QuantityValue<dimensions::MassDensity> central_density,
+        dimensions::QuantityValue<dimensions::GravitationalConstant> gravitational_constant) {
+        if (!finite_surface(solution) || !positive_finite(central_enthalpy.value()) || !positive_finite(central_density.value()) || !positive_finite(gravitational_constant.value())) {
+            return std::unexpected("Lane–Emden enthalpy scaling needs a finite surface and positive finite enthalpy, central density, and G.");
+        }
+        const long double alpha = std::sqrt(central_enthalpy.value() / (four_pi * gravitational_constant.value() * central_density.value()));
+        return make_scales(solution, alpha, central_density.value(), gravitational_constant.value());
+    }
+
+    std::expected<PhysicalScales, std::string> scale_from_polytropic_constant(
+        const Solution &solution, double polytropic_constant,
+        dimensions::QuantityValue<dimensions::MassDensity> central_density,
+        dimensions::QuantityValue<dimensions::GravitationalConstant> gravitational_constant) {
+        if (!finite_surface(solution) || solution.n <= 0.0 || !positive_finite(polytropic_constant) || !positive_finite(central_density.value()) ||
+            !positive_finite(gravitational_constant.value())) {
+            return std::unexpected("Lane–Emden EOS scaling needs n > 0, a finite surface, and positive finite K, central density, and G.");
+        }
+        const double enthalpy = static_cast<double>((solution.n + 1.0L) * polytropic_constant * std::pow(static_cast<long double>(central_density.value()), 1.0L / solution.n));
+        auto scales = scale_from_enthalpy(solution, dimensions::QuantityValue<dimensions::SpecificEnthalpy>{enthalpy}, central_density, gravitational_constant);
+        if (scales) {
+            scales->polytropic_constant = polytropic_constant;
+        }
+        return scales;
+    }
+
+    std::expected<PhysicalSample, std::string> evaluate_physical(
+        const Solution &solution, const PhysicalScales &scales, dimensions::QuantityValue<dimensions::Length> radius) {
+        if (!finite_surface(solution) || solution.n != scales.n || !std::isfinite(radius.value()) || radius.value() < 0.0 ||
+            !positive_finite(scales.radial_scale) || !positive_finite(scales.radius.value()) || !positive_finite(scales.mass.value()) ||
+            !positive_finite(scales.central_density.value()) || !positive_finite(scales.central_enthalpy.value()) || !positive_finite(scales.gravitational_constant.value())) {
+            return std::unexpected("Lane–Emden physical evaluation requires matching finite-surface scales and a finite nonnegative radius.");
+        }
+        const long double gravitational_mass = static_cast<long double>(scales.gravitational_constant.value()) * scales.mass.value();
+        double density = 0.0;
+        double pressure = 0.0;
+        double enthalpy = 0.0;
+        double mass = scales.mass.value();
+        double potential;
+        double gravity;
+        if (radius.value() >= scales.radius.value()) {
+            potential = static_cast<double>(-gravitational_mass / radius.value());
+            gravity = static_cast<double>(gravitational_mass / radius.value() / radius.value());
+        } else {
+            const double xi = std::min(radius.value() / scales.radial_scale, solution.surface->xi);
+            const auto sample = evaluate(solution, xi);
+            if (!sample) {
+                return std::unexpected(sample.error());
+            }
+            const double theta = std::max(0.0, sample->theta);
+            density = scales.central_density.value() * std::pow(theta, solution.n);
+            enthalpy = scales.central_enthalpy.value() * theta;
+            pressure = static_cast<double>(static_cast<long double>(density) * enthalpy / (solution.n + 1.0L));
+            potential = static_cast<double>(-gravitational_mass / scales.radius.value() - enthalpy);
+            gravity = -scales.central_enthalpy.value() / scales.radial_scale * sample->derivative;
+            mass = static_cast<double>(-four_pi * scales.central_density.value() * std::pow(static_cast<long double>(scales.radial_scale), 3) * xi * xi * sample->derivative);
+        }
+        if (!std::isfinite(density) || !std::isfinite(pressure) || !std::isfinite(enthalpy) || !std::isfinite(potential) || !std::isfinite(gravity) || !std::isfinite(mass)) {
+            return std::unexpected("Lane–Emden physical fields overflow in the supplied units.");
+        }
+        return PhysicalSample{
+            .density = dimensions::QuantityValue<dimensions::MassDensity>{density},
+            .pressure = dimensions::QuantityValue<dimensions::Pressure>{pressure},
+            .enthalpy = dimensions::QuantityValue<dimensions::SpecificEnthalpy>{enthalpy},
+            .potential = dimensions::QuantityValue<dimensions::GravitationalPotential>{potential},
+            .enclosed_mass = dimensions::QuantityValue<dimensions::Mass>{mass},
+            .gravity = dimensions::QuantityValue<dimensions::Acceleration>{gravity}
+        };
+    }
+}
diff --git a/src/lib/serif/utils/laneemden/solution.cpp b/src/lib/serif/utils/laneemden/solution.cpp
new file mode 100644
index 00000000..f5bf202b
--- /dev/null
+++ b/src/lib/serif/utils/laneemden/solution.cpp
@@ -0,0 +1,168 @@
+#include "serif/utils/laneemden/solution.hpp"
+
+#include <algorithm>
+#include <array>
+#include <cmath>
+#include <limits>
+
+namespace serif::utils::laneemden {
+    namespace {
+        /** @brief Evaluate Boudreaux's coupled a_k and c_k series iteratively, avoiding exponential recursive work. */
+        Sample central_sample(double n, double xi, std::size_t terms) {
+            std::array<double, 64> coefficients{};
+            std::array<double, 64> powers{};
+            coefficients[0] = 1.0;
+            powers[0] = 1.0;
+            for (std::size_t order = 1; order < terms; ++order) {
+                if (order >= 2) {
+                    coefficients[order] = -powers[order - 2] / static_cast<double>(order * (order + 1));
+                }
+                double sum = 0.0;
+                for (std::size_t k = 1; k <= order; ++k) {
+                    sum += (static_cast<double>(k) * (n + 1.0) - static_cast<double>(order)) * coefficients[k] * powers[order - k];
+                }
+                powers[order] = sum / static_cast<double>(order);
+            }
+
+            double theta = coefficients[terms - 1];
+            double derivative = 0.0;
+            for (std::size_t order = terms - 1; order > 0; --order) {
+                derivative = derivative * xi + theta;
+                theta = theta * xi + coefficients[order - 1];
+            }
+            return {.xi = xi, .theta = theta, .derivative = derivative};
+        }
+
+        /** @brief Evaluate the regular positive-density equation, using vacuum continuation only while locating the first zero. */
+        Sample slope(double n, const Sample &sample) {
+            const double source = std::pow(std::max(0.0, sample.theta), n);
+            return {.xi = 1.0, .theta = sample.derivative, .derivative = -2.0 * sample.derivative / sample.xi - source};
+        }
+
+        /** @brief Construct a Runge–Kutta stage with both coupled variables evaluated at the same radius. */
+        Sample stage(const Sample &initial, const Sample &rate, double step) {
+            return {
+                .xi = initial.xi + step,
+                .theta = initial.theta + step * rate.theta,
+                .derivative = initial.derivative + step * rate.derivative
+            };
+        }
+
+        /** @brief Advance the coupled Lane–Emden system with classical fourth-order Runge–Kutta. */
+        Sample advance(double n, const Sample &initial, double step) {
+            const Sample first = slope(n, initial);
+            const Sample second = slope(n, stage(initial, first, 0.5 * step));
+            const Sample third = slope(n, stage(initial, second, 0.5 * step));
+            const Sample fourth = slope(n, stage(initial, third, step));
+            return {
+                .xi = initial.xi + step,
+                .theta = initial.theta + step * (first.theta + 2.0 * second.theta + 2.0 * third.theta + fourth.theta) / 6.0,
+                .derivative = initial.derivative + step * (first.derivative + 2.0 * second.derivative + 2.0 * third.derivative + fourth.derivative) / 6.0
+            };
+        }
+
+        /** @brief Reject numerical states that have overflowed or left the regular positive branch. */
+        bool valid_interior(const Sample &sample) {
+            return std::isfinite(sample.xi) && std::isfinite(sample.theta) && std::isfinite(sample.derivative) &&
+                sample.theta > 0.0 && sample.derivative < 0.0;
+        }
+
+        /** @brief Resolve the first zero from a short bracket without raising negative theta to a fractional power. */
+        std::expected<Sample, std::string> locate_surface(double n, const Sample &initial, double width, double tolerance) {
+            double lower = 0.0;
+            double upper = width;
+            Sample result = advance(n, initial, upper);
+            if (!std::isfinite(result.theta) || result.theta > 0.0) {
+                return std::unexpected("Lane–Emden surface could not be bracketed; reduce the integration step.");
+            }
+            for (std::size_t iteration = 0; iteration < 80 && upper - lower > tolerance; ++iteration) {
+                const double middle = lower + 0.5 * (upper - lower);
+                result = advance(n, initial, middle);
+                if (result.theta > 0.0) {
+                    lower = middle;
+                } else {
+                    upper = middle;
+                }
+            }
+            result = advance(n, initial, lower + 0.5 * (upper - lower));
+            result.theta = 0.0;
+            if (!std::isfinite(result.derivative) || result.derivative >= 0.0) {
+                return std::unexpected("Lane–Emden surface derivative is not finite and negative; reduce the integration step.");
+            }
+            return result;
+        }
+    }
+
+    std::expected<Solution, std::string> solve(const double n, IntegrationOptions options) {
+        if (!std::isfinite(n) || n < 0.0) {
+            return std::unexpected("Lane–Emden index must be finite and nonnegative.");
+        }
+        if (!std::isfinite(options.step) || options.step <= 0.0 || !std::isfinite(options.maximum_radius) || options.maximum_radius <= 0.0 ||
+            !std::isfinite(options.surface_tolerance) || options.surface_tolerance <= 0.0 || options.surface_tolerance >= 0.01 ||
+            options.maximum_steps == 0 || options.central_series_terms < 6 || options.central_series_terms > 64) {
+            return std::unexpected("Lane–Emden integration needs positive finite step and maximum radius, a surface tolerance in (0, 0.01), a nonzero step budget, and 6 to 64 series terms.");
+        }
+
+        Solution solution{.n = n, .samples = {{.xi = 0.0, .theta = 1.0, .derivative = 0.0}}, .surface = std::nullopt};
+        const double first_radius = std::min({options.step, options.maximum_radius, 0.01 / std::sqrt(std::max(1.0, n))});
+        Sample current = central_sample(n, first_radius, options.central_series_terms);
+        if (!valid_interior(current)) {
+            return std::unexpected("Lane–Emden central series is not representable for the supplied index and step.");
+        }
+        solution.samples.push_back(current);
+
+        for (std::size_t iteration = 0; current.xi < options.maximum_radius; ++iteration) {
+            if (iteration >= options.maximum_steps) {
+                return std::unexpected("Lane–Emden integration exhausted maximum_steps before reaching the surface or maximum_radius.");
+            }
+            const double distance = -current.theta / current.derivative;
+            const double tolerance = options.surface_tolerance * std::max(1.0, current.xi);
+            if (distance <= tolerance && current.xi + 2.0 * distance <= options.maximum_radius) {
+                const double root_tolerance = std::max(4.0 * std::numeric_limits<double>::epsilon() * current.xi, tolerance * 0.01);
+                const auto surface = locate_surface(n, current, 2.0 * distance, root_tolerance);
+                if (!surface) {
+                    return std::unexpected(surface.error());
+                }
+                solution.surface = *surface;
+                solution.samples.push_back(*surface);
+                return solution;
+            }
+            // Shorten only the final approach; no accepted sample extends past the first zero.
+            const double step = std::min({options.step, 0.5 * distance, options.maximum_radius - current.xi});
+            if (current.xi + step == current.xi) {
+                return std::unexpected("Lane–Emden step cannot advance the radius at machine precision; relax surface_tolerance.");
+            }
+            current = advance(n, current, step);
+            if (!valid_interior(current)) {
+                return std::unexpected("Lane–Emden integration left the finite positive branch; reduce the integration step.");
+            }
+            solution.samples.push_back(current);
+        }
+        return solution;
+    }
+
+    std::expected<Sample, std::string> evaluate(const Solution &solution, double xi) {
+        if (!std::isfinite(xi) || xi < 0.0 || solution.samples.empty() || xi > solution.samples.back().xi) {
+            return std::unexpected("Lane–Emden evaluation radius must lie within the stored finite dimensionless profile.");
+        }
+        const auto upper = std::lower_bound(solution.samples.begin(), solution.samples.end(), xi,
+            [](const Sample &sample, double radius) { return sample.xi < radius; });
+        if (upper == solution.samples.begin() || upper->xi == xi) {
+            return *upper;
+        }
+        const Sample &left = *(upper - 1);
+        const Sample &right = *upper;
+        const double width = right.xi - left.xi;
+        if (!std::isfinite(width) || width <= 0.0) {
+            return std::unexpected("Lane–Emden evaluation requires a profile with increasing finite radii produced by solve.");
+        }
+        const double position = (xi - left.xi) / width;
+        const double quadratic = 3.0 * (right.theta - left.theta) - width * (2.0 * left.derivative + right.derivative);
+        const double cubic = 2.0 * (left.theta - right.theta) + width * (left.derivative + right.derivative);
+        return Sample{
+            .xi = xi,
+            .theta = left.theta + position * (width * left.derivative + position * (quadratic + position * cubic)),
+            .derivative = left.derivative + position * (2.0 * quadratic + 3.0 * position * cubic) / width
+        };
+    }
+}