diff --git a/include/NeoFOAM/dsl/coeff.hpp b/include/NeoFOAM/dsl/coeff.hpp
new file mode 100644
index 000000000..cd91ecdb7
--- /dev/null
+++ b/include/NeoFOAM/dsl/coeff.hpp
@@ -0,0 +1,110 @@
+// SPDX-License-Identifier: MIT
+// SPDX-FileCopyrightText: 2023 NeoFOAM authors
+#pragma once
+
+namespace NeoFOAM::dsl
+{
+
+/**
+ * @class Coeff
+ * @brief A class that represents a coefficient for the NeoFOAM dsl.
+ *
+ * This class stores a single scalar coefficient and optionally span of values.
+ * It is used to delay the evaluation of a scalar multiplication with a field to
+ * avoid the creation of a temporary field copy.
+ * It provides an indexing operator `operator[]` that returns the evaluated value at the specified
+ * index.
+ */
+class Coeff
+{
+
+public:
+
+    Coeff() : coeff_(1.0), span_(), hasSpan_(false) {}
+
+    Coeff(scalar value) : coeff_(value), span_(), hasSpan_(false) {}
+
+    Coeff(scalar coeff, const Field<scalar>& field)
+        : coeff_(coeff), span_(field.span()), hasSpan_(true)
+    {}
+
+    Coeff(const Field<scalar>& field) : coeff_(1.0), span_(field.span()), hasSpan_(true) {}
+
+    KOKKOS_INLINE_FUNCTION
+    scalar operator[](const size_t i) const { return (hasSpan_) ? span_[i] * coeff_ : coeff_; }
+
+    bool hasSpan() { return hasSpan_; }
+
+    std::span<const scalar> span() { return span_; }
+
+    Coeff& operator*=(scalar rhs)
+    {
+        coeff_ *= rhs;
+        return *this;
+    }
+
+    Coeff& operator*=(const Coeff& rhs)
+    {
+        if (hasSpan_ && rhs.hasSpan_)
+        {
+            NF_ERROR_EXIT("Not implemented");
+        }
+
+        if (!hasSpan_ && rhs.hasSpan_)
+        {
+            // Take over the span
+            span_ = rhs.span_;
+            hasSpan_ = true;
+        }
+
+        return this->operator*=(rhs.coeff_);
+    }
+
+
+private:
+
+    scalar coeff_;
+
+    std::span<const scalar> span_;
+
+    bool hasSpan_;
+};
+
+
+namespace detail
+
+{
+/* @brief function to force evaluation to a field, the field will be resized to hold either a
+ * single value or the full field
+ *
+ * @param field to store the result
+ */
+void toField(Coeff& coeff, Field<scalar>& rhs)
+{
+    if (coeff.hasSpan())
+    {
+        rhs.resize(coeff.span().size());
+        fill(rhs, 1.0);
+        auto rhsSpan = rhs.span();
+        // otherwise we are unable to capture values in the lambda
+        parallelFor(
+            rhs.exec(), rhs.range(), KOKKOS_LAMBDA(const size_t i) { rhsSpan[i] *= coeff[i]; }
+        );
+    }
+    else
+    {
+        rhs.resize(1);
+        fill(rhs, coeff[0]);
+    }
+}
+
+}
+
+inline Coeff operator*(const Coeff& lhs, const Coeff& rhs)
+{
+    Coeff result = lhs;
+    result *= rhs;
+    return result;
+}
+
+} // namespace NeoFOAM::dsl
diff --git a/include/NeoFOAM/dsl/ddt.hpp b/include/NeoFOAM/dsl/ddt.hpp
new file mode 100644
index 000000000..06a4b5ba3
--- /dev/null
+++ b/include/NeoFOAM/dsl/ddt.hpp
@@ -0,0 +1,50 @@
+// SPDX-License-Identifier: MIT
+//
+// SPDX-FileCopyrightText: 2023 NeoFOAM authors
+
+#pragma once
+
+#include "NeoFOAM/fields/field.hpp"
+#include "NeoFOAM/finiteVolume/cellCentred.hpp"
+#include "NeoFOAM/finiteVolume/cellCentred/fields/volumeField.hpp"
+
+namespace NeoFOAM::dsl::temporal
+{
+
+
+// TODO add free factory function
+template<typename FieldType>
+class Ddt : public OperatorMixin<FieldType>
+{
+
+public:
+
+    Ddt(FieldType& field) : OperatorMixin<FieldType>(field.exec(), field, Operator::Type::Temporal)
+    {}
+
+    std::string getName() const { return "TimeOperator"; }
+
+    void explicitOperation([[maybe_unused]] Field<scalar>& source, [[maybe_unused]] scalar scale)
+    {
+        NF_ERROR_EXIT("Not implemented");
+    }
+
+    void implicitOperation([[maybe_unused]] Field<scalar>& phi)
+    {
+        NF_ERROR_EXIT("Not implemented");
+    }
+
+private:
+};
+
+// see
+// https://github.com/exasim-project/NeoFOAM/blob/dsl/operatorIntergration/include/NeoFOAM/finiteVolume/cellCentred/operators/explicitOperators/expOp.hpp
+
+template<typename FieldType>
+Ddt<FieldType> ddt(FieldType& in)
+{
+    return Ddt(in);
+};
+
+
+} // namespace NeoFOAM
diff --git a/include/NeoFOAM/dsl/expression.hpp b/include/NeoFOAM/dsl/expression.hpp
new file mode 100644
index 000000000..90808f1d2
--- /dev/null
+++ b/include/NeoFOAM/dsl/expression.hpp
@@ -0,0 +1,170 @@
+// SPDX-License-Identifier: MIT
+// SPDX-FileCopyrightText: 2023-2024 NeoFOAM authors
+#pragma once
+
+#include <iostream>
+#include <memory>
+#include <vector>
+#include <utility>
+
+#include "NeoFOAM/core/primitives/scalar.hpp"
+#include "NeoFOAM/fields/field.hpp"
+#include "NeoFOAM/dsl/operator.hpp"
+#include "NeoFOAM/core/error.hpp"
+
+#include "NeoFOAM/finiteVolume/cellCentred/fields/volumeField.hpp"
+
+namespace NeoFOAM::dsl
+{
+
+
+class Expression
+{
+public:
+
+    Expression(const Executor& exec)
+        : exec_(exec), temporalOperators_(), implicitOperators_(), explicitOperators_()
+    {}
+
+    /* @brief perform all explicit operation and accumulate the result */
+    Field<scalar> explicitOperation(size_t nCells) const
+    {
+        Field<scalar> source(exec_, nCells, 0.0);
+        return explicitOperation(source);
+    }
+
+    /* @brief perform all explicit operation and accumulate the result */
+    Field<scalar> explicitOperation(Field<scalar>& source) const
+    {
+        for (auto& Operator : explicitOperators_)
+        {
+            Operator.explicitOperation(source);
+        }
+        return source;
+    }
+
+    void addOperator(const Operator& Operator)
+    {
+        switch (Operator.getType())
+        {
+        case Operator::Type::Temporal:
+            temporalOperators_.push_back(Operator);
+            break;
+        case Operator::Type::Implicit:
+            implicitOperators_.push_back(Operator);
+            break;
+        case Operator::Type::Explicit:
+            explicitOperators_.push_back(Operator);
+            break;
+        }
+    }
+
+    void addExpression(const Expression& equation)
+    {
+        for (auto& Operator : equation.temporalOperators_)
+        {
+            temporalOperators_.push_back(Operator);
+        }
+        for (auto& Operator : equation.implicitOperators_)
+        {
+            implicitOperators_.push_back(Operator);
+        }
+        for (auto& Operator : equation.explicitOperators_)
+        {
+            explicitOperators_.push_back(Operator);
+        }
+    }
+
+
+    /* @brief getter for the total number of terms in the equation */
+    size_t size() const
+    {
+        return temporalOperators_.size() + implicitOperators_.size() + explicitOperators_.size();
+    }
+
+    // getters
+    const std::vector<Operator>& temporalOperators() const { return temporalOperators_; }
+
+    const std::vector<Operator>& implicitOperators() const { return implicitOperators_; }
+
+    const std::vector<Operator>& explicitOperators() const { return explicitOperators_; }
+
+    std::vector<Operator>& temporalOperators() { return temporalOperators_; }
+
+    std::vector<Operator>& implicitOperators() { return implicitOperators_; }
+
+    std::vector<Operator>& explicitOperators() { return explicitOperators_; }
+
+    const Executor& exec() const { return exec_; }
+
+private:
+
+    const Executor exec_;
+
+    std::vector<Operator> temporalOperators_;
+
+    std::vector<Operator> implicitOperators_;
+
+    std::vector<Operator> explicitOperators_;
+};
+
+Expression operator+(Expression lhs, const Expression& rhs)
+{
+    lhs.addExpression(rhs);
+    return lhs;
+}
+
+Expression operator+(Expression lhs, const Operator& rhs)
+{
+    lhs.addOperator(rhs);
+    return lhs;
+}
+
+Expression operator+(const Operator& lhs, const Operator& rhs)
+{
+    Expression expr(lhs.exec());
+    expr.addOperator(lhs);
+    expr.addOperator(rhs);
+    return expr;
+}
+
+Expression operator*(scalar scale, const Expression& es)
+{
+    Expression expr(es.exec());
+    for (const auto& Operator : es.temporalOperators())
+    {
+        expr.addOperator(scale * Operator);
+    }
+    for (const auto& Operator : es.implicitOperators())
+    {
+        expr.addOperator(scale * Operator);
+    }
+    for (const auto& Operator : es.explicitOperators())
+    {
+        expr.addOperator(scale * Operator);
+    }
+    return expr;
+}
+
+Expression operator-(Expression lhs, const Expression& rhs)
+{
+    lhs.addExpression(-1.0 * rhs);
+    return lhs;
+}
+
+Expression operator-(Expression lhs, const Operator& rhs)
+{
+    lhs.addOperator(-1.0 * rhs);
+    return lhs;
+}
+
+Expression operator-(const Operator& lhs, const Operator& rhs)
+{
+    Expression expr(lhs.exec());
+    expr.addOperator(lhs);
+    expr.addOperator(Coeff(-1) * rhs);
+    return expr;
+}
+
+
+} // namespace NeoFOAM::dsl
diff --git a/include/NeoFOAM/dsl/operator.hpp b/include/NeoFOAM/dsl/operator.hpp
new file mode 100644
index 000000000..345432fa9
--- /dev/null
+++ b/include/NeoFOAM/dsl/operator.hpp
@@ -0,0 +1,259 @@
+// SPDX-License-Identifier: MIT
+// SPDX-FileCopyrightText: 2023-2024 NeoFOAM authors
+#pragma once
+
+#include <iostream>
+#include <memory>
+#include <type_traits>
+#include <utility>
+#include <concepts>
+
+#include "NeoFOAM/core/primitives/scalar.hpp"
+#include "NeoFOAM/fields/field.hpp"
+#include "NeoFOAM/core/input.hpp"
+#include "NeoFOAM/dsl/coeff.hpp"
+#include "NeoFOAM/finiteVolume/cellCentred/fields/volumeField.hpp"
+
+namespace NeoFOAM::dsl
+{
+
+template<typename T>
+concept HasTemporalOperator = requires(T t) {
+    {
+        t.temporalOperation(std::declval<Field<scalar>&>(), std::declval<scalar>())
+    } -> std::same_as<void>; // Adjust return type and arguments as needed
+};
+
+template<typename T>
+concept HasExplicitOperator = requires(T t) {
+    {
+        t.explicitOperation(std::declval<Field<scalar>&>())
+    } -> std::same_as<void>; // Adjust return type and arguments as needed
+};
+
+/* @class Operator
+ * @brief A class to represent a operator in NeoFOAMs dsl
+ *
+ * The design here is based on the type erasure design pattern
+ * see https://www.youtube.com/watch?v=4eeESJQk-mw
+ *
+ * Motivation for using type erasure is that concrete implementation
+ * of Operators e.g Divergence, Laplacian, etc can be stored in a vector of
+ * Operators
+ *
+ * @ingroup dsl
+ */
+class Operator
+{
+public:
+
+    enum class Type
+    {
+        Temporal,
+        Implicit,
+        Explicit
+    };
+
+    template<typename T>
+    Operator(T cls) : model_(std::make_unique<OperatorModel<T>>(std::move(cls)))
+    {}
+
+    Operator(const Operator& eqnOperator) : model_ {eqnOperator.model_->clone()} {}
+
+    Operator(Operator&& eqnOperator) : model_ {std::move(eqnOperator.model_)} {}
+
+    void explicitOperation(Field<scalar>& source) const { model_->explicitOperation(source); }
+
+    void temporalOperation(Field<scalar>& field) { model_->temporalOperation(field); }
+
+    /* returns the fundamental type of an operator, ie explicit, implicit, temporal */
+    Operator::Type getType() const { return model_->getType(); }
+
+    Coeff& getCoefficient() { return model_->getCoefficient(); }
+
+    Coeff getCoefficient() const { return model_->getCoefficient(); }
+
+    /* get the corresponding field size over which the operator operates */
+    // size_t getSize() const { return model_->getSize(); }
+
+    /* @brief Given an input this function reads required coeffs */
+    void build(const Input& input) { model_->build(input); }
+
+    /* @brief Get the executor */
+    const Executor& exec() const { return model_->exec(); }
+
+
+private:
+
+    /* @brief Base class defining the concept of a term. This effectively
+     * defines what functions need to be implemented by a concrete Operator implementation
+     * */
+    struct OperatorConcept
+    {
+        virtual ~OperatorConcept() = default;
+
+        virtual void explicitOperation(Field<scalar>& source) const = 0;
+
+        virtual void temporalOperation(Field<scalar>& field) = 0;
+
+        /* @brief Given an input this function reads required coeffs */
+        virtual void build(const Input& input) = 0;
+
+        /* returns the name of the operator */
+        virtual std::string getName() const = 0;
+
+        /* returns the fundamental type of an operator, ie explicit, implicit, temporal */
+        virtual Operator::Type getType() const = 0;
+
+        // virtual size_t getSize() const = 0;
+
+        /* @brief get the associated coefficient for this term */
+        virtual Coeff& getCoefficient() = 0;
+
+        /* @brief get the associated coefficient for this term */
+        virtual Coeff getCoefficient() const = 0;
+
+        /* @brief Get the executor */
+        virtual const Executor& exec() const = 0;
+
+        // The Prototype Design Pattern
+        virtual std::unique_ptr<OperatorConcept> clone() const = 0;
+    };
+
+    // Templated derived class to implement the type-specific behavior
+    template<typename ConcreteOperatorType>
+    struct OperatorModel : OperatorConcept
+    {
+        /* @brief build with concrete operator */
+        OperatorModel(ConcreteOperatorType concreteOp) : concreteOp_(std::move(concreteOp)) {}
+
+        /* returns the name of the operator */
+        std::string getName() const override { return concreteOp_.getName(); }
+
+        virtual void explicitOperation(Field<scalar>& source) const override
+        {
+            if constexpr (HasExplicitOperator<ConcreteOperatorType>)
+            {
+                concreteOp_.explicitOperation(source);
+            }
+        }
+
+        virtual void temporalOperation(Field<scalar>& field) override
+        {
+            if constexpr (HasTemporalOperator<ConcreteOperatorType>)
+            {
+                concreteOp_.temporalOperation(field);
+            }
+        }
+
+        /* @brief Given an input this function reads required coeffs */
+        virtual void build(const Input& input) override { concreteOp_.build(input); }
+
+        /* returns the fundamental type of an operator, ie explicit, implicit, temporal */
+        Operator::Type getType() const override { return concreteOp_.getType(); }
+
+        /* @brief Get the executor */
+        const Executor& exec() const override { return concreteOp_.exec(); }
+
+        /* @brief get the associated coefficient for this term */
+        virtual Coeff& getCoefficient() override { return concreteOp_.getCoefficient(); }
+
+        /* @brief get the associated coefficient for this term */
+        virtual Coeff getCoefficient() const override { return concreteOp_.getCoefficient(); }
+
+        /* @brief get the associated coefficient for this term */
+        // virtual size_t getSize() const override { return concreteOp_.getSize(); }
+
+        // The Prototype Design Pattern
+        std::unique_ptr<OperatorConcept> clone() const override
+        {
+            return std::make_unique<OperatorModel>(*this);
+        }
+
+        ConcreteOperatorType concreteOp_;
+    };
+
+    std::unique_ptr<OperatorConcept> model_;
+};
+
+
+auto operator*(scalar scalarCoeff, const Operator& rhs)
+{
+    Operator result = rhs;
+    result.getCoefficient() *= scalarCoeff;
+    return result;
+}
+
+auto operator*(const Field<scalar>& coeffField, const Operator& rhs)
+{
+    Operator result = rhs;
+    result.getCoefficient() *= Coeff(coeffField);
+    return result;
+}
+
+auto operator*(const Coeff& coeff, const Operator& rhs)
+{
+    Operator result = rhs;
+    result.getCoefficient() *= coeff;
+    return result;
+}
+
+template<typename CoeffFunction>
+    requires std::invocable<CoeffFunction&, size_t>
+Operator operator*([[maybe_unused]] CoeffFunction coeffFunc, const Operator& lhs)
+{
+    // TODO implement
+    NF_ERROR_EXIT("Not implemented");
+    Operator result = lhs;
+    // if (!result.getCoefficient().useSpan)
+    // {
+    //     result.setField(std::make_shared<Field<scalar>>(result.exec(), result.nCells(), 1.0));
+    // }
+    // map(result.exec(), result.getCoefficient().values, scaleFunc);
+    return result;
+}
+
+/* @class OperatorMixin
+ * @brief A mixin class to represent simplify implementations of concrete operators
+ * in NeoFOAMs dsl
+ *
+ *
+ * @ingroup dsl
+ */
+template<typename FieldType>
+class OperatorMixin
+{
+
+public:
+
+    OperatorMixin(const Executor exec, FieldType& field, Operator::Type type)
+        : exec_(exec), coeffs_(), field_(field), type_(type) {};
+
+    Operator::Type getType() const { return type_; }
+
+    virtual ~OperatorMixin() = default;
+
+    virtual const Executor& exec() const final { return exec_; }
+
+    Coeff& getCoefficient() { return coeffs_; }
+
+    const Coeff& getCoefficient() const { return coeffs_; }
+
+    FieldType& getField() { return field_; }
+
+    const FieldType& getField() const { return field_; }
+
+    /* @brief Given an input this function reads required coeffs */
+    void build([[maybe_unused]] const Input& input) {}
+
+protected:
+
+    const Executor exec_; //!< Executor associated with the field. (CPU, GPU, openMP, etc.)
+
+    Coeff coeffs_;
+
+    FieldType& field_;
+
+    Operator::Type type_;
+};
+} // namespace NeoFOAM::dsl
diff --git a/include/NeoFOAM/dsl/solver.hpp b/include/NeoFOAM/dsl/solver.hpp
new file mode 100644
index 000000000..2f84fe57b
--- /dev/null
+++ b/include/NeoFOAM/dsl/solver.hpp
@@ -0,0 +1,55 @@
+// SPDX-License-Identifier: MIT
+// SPDX-FileCopyrightText: 2023-2024 NeoFOAM authors
+#pragma once
+
+#include <iostream>
+#include <memory>
+#include <type_traits>
+#include <utility>
+#include <concepts>
+
+#include "NeoFOAM/core/primitives/scalar.hpp"
+#include "NeoFOAM/fields/field.hpp"
+#include "NeoFOAM/core/input.hpp"
+#include "NeoFOAM/dsl/expression.hpp"
+#include "NeoFOAM/timeIntegration/timeIntegration.hpp"
+
+
+namespace NeoFOAM::dsl
+{
+
+/* @brief solve an expression
+ *
+ * @param exp - Expression which is to be solved/updated.
+ * @param solution - Solution field, where the solution will be 'written to'.
+ * @param dt - time step for the temporal integration
+ * @param fvSchemes - Dictionary containing spatial operator and time  integration properties
+ * @param fvSolution - Dictionary containing linear solver properties
+ */
+template<typename FieldType>
+void solve(
+    Expression& exp,
+    FieldType& solution,
+    scalar dt,
+    [[maybe_unused]] const Dictionary& fvSchemes,
+    [[maybe_unused]] const Dictionary& fvSolution
+)
+{
+    if (exp.temporalOperators().size() == 0 && exp.implicitOperators().size() == 0)
+    {
+        NF_ERROR_EXIT("No temporal or implicit terms to solve.");
+    }
+    if (exp.temporalOperators().size() > 0)
+    {
+        // integrate equations in time
+        TimeIntegration<FieldType> timeIntegrator(fvSchemes.subDict("ddtSchemes"));
+        timeIntegrator.solve(exp, solution, dt);
+    }
+    else
+    {
+        NF_ERROR_EXIT("Not implemented.");
+        // solve sparse matrix system
+    }
+}
+
+} // namespace NeoFOAM::dsl
diff --git a/include/NeoFOAM/fields/field.hpp b/include/NeoFOAM/fields/field.hpp
index c5e02ab51..fd5ffd691 100644
--- a/include/NeoFOAM/fields/field.hpp
+++ b/include/NeoFOAM/fields/field.hpp
@@ -79,7 +79,7 @@ class Field
     {
         void* ptr = nullptr;
         std::visit(
-            [this, &ptr, size](const auto& concreteExec)
+            [&ptr, size](const auto& concreteExec)
             { ptr = concreteExec.alloc(size * sizeof(ValueType)); },
             exec_
         );
@@ -413,7 +413,13 @@ class Field
     void validateOtherField(const Field<ValueType>& rhs) const
     {
         NF_DEBUG_ASSERT(size() == rhs.size(), "Fields are not the same size.");
-        NF_DEBUG_ASSERT(exec() == rhs.exec(), "Executors are not the same.");
+
+        std::string execName = std::visit([](auto e) { return e.print(); }, exec());
+        std::string rhsExecName = std::visit([](auto e) { return e.print(); }, rhs.exec());
+        NF_DEBUG_ASSERT(
+            exec() == rhs.exec(),
+            "Executors: " + execName + " and " + rhsExecName + " are not the same."
+        );
     }
 };
 
diff --git a/include/NeoFOAM/finiteVolume/cellCentred/fields/geometricField.hpp b/include/NeoFOAM/finiteVolume/cellCentred/fields/geometricField.hpp
index c02ee804a..d2fb1dea4 100644
--- a/include/NeoFOAM/finiteVolume/cellCentred/fields/geometricField.hpp
+++ b/include/NeoFOAM/finiteVolume/cellCentred/fields/geometricField.hpp
@@ -78,6 +78,13 @@ class GeometricFieldMixin
      */
     Field<ValueType>& internalField() { return field_.internalField(); }
 
+    /**
+     * @brief Returns the size of the internal field
+     *
+     * @return The size of the internal field
+     */
+    size_t size() const { return field_.internalField().size(); }
+
     /**
      * @brief Returns a const reference to the boundary field.
      *
diff --git a/include/NeoFOAM/timeIntegration/forwardEuler.hpp b/include/NeoFOAM/timeIntegration/forwardEuler.hpp
new file mode 100644
index 000000000..0dda395a1
--- /dev/null
+++ b/include/NeoFOAM/timeIntegration/forwardEuler.hpp
@@ -0,0 +1,56 @@
+// SPDX-License-Identifier: MIT
+//
+// SPDX-FileCopyrightText: 2023 NeoFOAM authors
+
+#pragma once
+
+#include <functional>
+
+#include "NeoFOAM/fields/field.hpp"
+#include "NeoFOAM/timeIntegration/timeIntegration.hpp"
+
+namespace NeoFOAM::dsl
+{
+
+template<typename SolutionType>
+class ForwardEuler :
+    public TimeIntegratorBase<SolutionType>::template Register<ForwardEuler<SolutionType>>
+{
+
+public:
+
+    using Base = TimeIntegratorBase<SolutionType>::template Register<ForwardEuler<SolutionType>>;
+
+    ForwardEuler(const Dictionary& dict) : Base(dict) {}
+
+    static std::string name() { return "forwardEuler"; }
+
+    static std::string doc() { return "first order time integration method"; }
+
+    static std::string schema() { return "none"; }
+
+    void solve(Expression& eqn, SolutionType& sol, scalar dt) const override
+    {
+        auto source = eqn.explicitOperation(sol.size());
+
+        sol.internalField() -= source * dt;
+        sol.correctBoundaryConditions();
+
+        // check if executor is GPU
+        if (std::holds_alternative<NeoFOAM::GPUExecutor>(eqn.exec()))
+        {
+            Kokkos::fence();
+        }
+    };
+
+    std::unique_ptr<TimeIntegratorBase<SolutionType>> clone() const override
+    {
+        return std::make_unique<ForwardEuler>(*this);
+    }
+};
+
+// unfortunately needs explicit instantiation
+template class ForwardEuler<finiteVolume::cellCentred::VolumeField<scalar>>;
+
+
+} // namespace NeoFOAM
diff --git a/include/NeoFOAM/timeIntegration/timeIntegration.hpp b/include/NeoFOAM/timeIntegration/timeIntegration.hpp
new file mode 100644
index 000000000..89b1feff9
--- /dev/null
+++ b/include/NeoFOAM/timeIntegration/timeIntegration.hpp
@@ -0,0 +1,76 @@
+// SPDX-License-Identifier: MIT
+//
+// SPDX-FileCopyrightText: 2023 NeoFOAM authors
+
+#pragma once
+
+#include <functional>
+
+#include "NeoFOAM/fields/field.hpp"
+#include "NeoFOAM/finiteVolume/cellCentred.hpp"
+#include "NeoFOAM/dsl/expression.hpp"
+
+namespace NeoFOAM::dsl
+{
+
+/* @class Factory class to create time integration method by a given name
+ * using NeoFOAMs runTimeFactory mechanism
+ */
+template<typename SolutionType>
+class TimeIntegratorBase :
+    public RuntimeSelectionFactory<TimeIntegratorBase<SolutionType>, Parameters<const Dictionary&>>
+{
+
+public:
+
+    static std::string name() { return "timeIntegrationFactory"; }
+
+    TimeIntegratorBase(const Dictionary& dict) : dict_(dict) {}
+
+    virtual ~TimeIntegratorBase() {}
+
+    virtual void solve(Expression& eqn, SolutionType& sol, scalar dt)
+        const = 0; // Pure virtual function for solving
+
+    // Pure virtual function for cloning
+    virtual std::unique_ptr<TimeIntegratorBase> clone() const = 0;
+
+protected:
+
+    const Dictionary& dict_;
+};
+
+/* @class Factory class to create time integration method by a given name
+ * using NeoFOAMs runTimeFactory mechanism
+ *
+ * @tparam SolutionType Type of the solution field eg, volumeField or just a plain Field
+ */
+template<typename SolutionType>
+class TimeIntegration
+{
+
+public:
+
+    TimeIntegration(const TimeIntegration& timeIntegrator)
+        : timeIntegratorStrategy_(timeIntegrator.timeIntegratorStrategy_->clone()) {};
+
+    TimeIntegration(TimeIntegration&& timeIntegrator)
+        : timeIntegratorStrategy_(std::move(timeIntegrator.timeIntegratorStrategy_)) {};
+
+    TimeIntegration(const Dictionary& dict)
+        : timeIntegratorStrategy_(
+            TimeIntegratorBase<SolutionType>::create(dict.get<std::string>("type"), dict)
+        ) {};
+
+    void solve(Expression& eqn, SolutionType& sol, scalar dt)
+    {
+        timeIntegratorStrategy_->solve(eqn, sol, dt);
+    }
+
+private:
+
+    std::unique_ptr<TimeIntegratorBase<SolutionType>> timeIntegratorStrategy_;
+};
+
+
+} // namespace NeoFOAM::dsl
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 80d8ffdce..a0dbcffef 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -6,7 +6,8 @@ list(APPEND CMAKE_MODULE_PATH ${Catch2_SOURCE_DIR}/extras)
 # has to be first since it adds the main target
 add_subdirectory(catch2)
 
-# This function creates unit tests. It allows the following keywords:
+# include "NeoFOAM/core/parallelAlgorithms.hpp" This function creates unit tests. It provides the
+# following keywords:
 #
 # * MPI_SIZE: the number of MPI processors to be used, defaults to 1 if not set
 # * COMMAND: the test command (same behavior as for CMake's add_test), defaults to the test name
diff --git a/test/dsl/CMakeLists.txt b/test/dsl/CMakeLists.txt
index 4935deff8..419310dd5 100644
--- a/test/dsl/CMakeLists.txt
+++ b/test/dsl/CMakeLists.txt
@@ -3,4 +3,6 @@
 
 neofoam_unit_test(coeff)
 neofoam_unit_test(operator)
-neofoam_unit_test(equation)
+neofoam_unit_test(expression)
+
+add_subdirectory(timeIntegration)
\ No newline at end of file
diff --git a/test/dsl/coeff.cpp b/test/dsl/coeff.cpp
index e7c33958e..f945a00ae 100644
--- a/test/dsl/coeff.cpp
+++ b/test/dsl/coeff.cpp
@@ -8,12 +8,11 @@
 #include <catch2/benchmark/catch_benchmark.hpp>
 
 #include "NeoFOAM/fields/field.hpp"
-#include "NeoFOAM/DSL/coeff.hpp"
+#include "NeoFOAM/dsl/coeff.hpp"
 
 using Field = NeoFOAM::Field<NeoFOAM::scalar>;
-using Coeff = NeoFOAM::DSL::Coeff;
-
-namespace detail = NeoFOAM::DSL::detail;
+using Coeff = NeoFOAM::dsl::Coeff;
+using namespace NeoFOAM::dsl;
 
 
 TEST_CASE("Coeff")
@@ -62,8 +61,8 @@ TEST_CASE("Coeff")
     {
         size_t size = 3;
 
-        NeoFOAM::Field<NeoFOAM::scalar> fieldA(exec, size, 0.0);
-        NeoFOAM::Field<NeoFOAM::scalar> fieldB(exec, size, 1.0);
+        Field fieldA(exec, size, 0.0);
+        Field fieldB(exec, size, 1.0);
 
         SECTION("span")
         {
@@ -82,8 +81,6 @@ TEST_CASE("Coeff")
 
         SECTION("scalar")
         {
-            NeoFOAM::Field<NeoFOAM::scalar> fieldA(exec, size, 0.0);
-
             Coeff coeff = Coeff(2.0);
             {
                 NeoFOAM::parallelFor(
@@ -99,8 +96,6 @@ TEST_CASE("Coeff")
 
         SECTION("span and scalar")
         {
-            NeoFOAM::Field<NeoFOAM::scalar> fieldA(exec, size, 0.0);
-            NeoFOAM::Field<NeoFOAM::scalar> fieldB(exec, size, 1.0);
             Coeff coeff {-5.0, fieldB};
             {
                 NeoFOAM::parallelFor(
diff --git a/test/dsl/common.hpp b/test/dsl/common.hpp
index 434c78e7f..7a7115905 100644
--- a/test/dsl/common.hpp
+++ b/test/dsl/common.hpp
@@ -8,17 +8,19 @@
 #include <catch2/benchmark/catch_benchmark.hpp>
 
 #include "NeoFOAM/fields/field.hpp"
-#include "NeoFOAM/DSL/coeff.hpp"
-#include "NeoFOAM/DSL/operator.hpp"
+#include "NeoFOAM/finiteVolume/cellCentred/fields/volumeField.hpp"
+#include "NeoFOAM/dsl/coeff.hpp"
+#include "NeoFOAM/dsl/operator.hpp"
+
+namespace fvcc = NeoFOAM::finiteVolume::cellCentred;
 
 using Field = NeoFOAM::Field<NeoFOAM::scalar>;
-using Coeff = NeoFOAM::DSL::Coeff;
-using Operator = NeoFOAM::DSL::Operator;
-using OperatorMixin = NeoFOAM::DSL::OperatorMixin;
+using Coeff = NeoFOAM::dsl::Coeff;
+using Operator = NeoFOAM::dsl::Operator;
 using Executor = NeoFOAM::Executor;
 using VolumeField = fvcc::VolumeField<NeoFOAM::scalar>;
+using OperatorMixin = NeoFOAM::dsl::OperatorMixin<VolumeField>;
 using BoundaryFields = NeoFOAM::BoundaryFields<NeoFOAM::scalar>;
-namespace fvcc = NeoFOAM::finiteVolume::cellCentred;
 
 /* A dummy implementation of a Operator
  * following the Operator interface */
@@ -27,7 +29,7 @@ class Dummy : public OperatorMixin
 
 public:
 
-    Dummy(const Executor& exec, VolumeField& field) : OperatorMixin(exec), field_(field) {}
+    Dummy(VolumeField& field) : OperatorMixin(field.exec(), field, Operator::Type::Explicit) {}
 
     void explicitOperation(Field& source) const
     {
@@ -42,18 +44,6 @@ class Dummy : public OperatorMixin
     }
 
     std::string getName() const { return "Dummy"; }
-
-    const VolumeField& volumeField() const { return field_; }
-
-    VolumeField& volumeField() { return field_; }
-
-    Operator::Type getType() const { return Operator::Type::Explicit; }
-
-    size_t getSize() const { return field_.internalField().size(); }
-
-private:
-
-    VolumeField& field_;
 };
 
 NeoFOAM::scalar getField(const NeoFOAM::Field<NeoFOAM::scalar>& source)
diff --git a/test/dsl/expression.cpp b/test/dsl/expression.cpp
new file mode 100644
index 000000000..c1d509e36
--- /dev/null
+++ b/test/dsl/expression.cpp
@@ -0,0 +1,52 @@
+// SPDX-License-Identifier: MIT
+// SPDX-FileCopyrightText: 2023-2024 NeoFOAM authors
+#define CATCH_CONFIG_RUNNER // Define this before including catch.hpp to create
+                            // a custom main
+#include "common.hpp"
+#include "NeoFOAM/dsl/expression.hpp"
+
+using Expression = NeoFOAM::dsl::Expression;
+
+TEST_CASE("Expression")
+{
+    NeoFOAM::Executor exec = GENERATE(
+        NeoFOAM::Executor(NeoFOAM::SerialExecutor {}),
+        NeoFOAM::Executor(NeoFOAM::CPUExecutor {}),
+        NeoFOAM::Executor(NeoFOAM::GPUExecutor {})
+    );
+
+    std::string execName = std::visit([](auto e) { return e.print(); }, exec);
+    auto mesh = NeoFOAM::createSingleCellMesh(exec);
+
+    const size_t size {1};
+    Field fA(exec, size, 2.0);
+    BoundaryFields bf(exec, mesh.nBoundaryFaces(), mesh.nBoundaries());
+
+    std::vector<fvcc::VolumeBoundary<NeoFOAM::scalar>> bcs {};
+    auto vf = VolumeField(exec, mesh, fA, bf, bcs);
+    auto fB = Field(exec, 1, 4.0);
+
+    auto a = Dummy(vf);
+    auto b = Dummy(vf);
+
+    SECTION("Create equation and perform explicitOperation on " + execName)
+    {
+        auto eqnA = a + b;
+        auto eqnB = fB * Dummy(vf) + 2 * Dummy(vf);
+        auto eqnC = Expression(2 * a - b);
+        auto eqnD = 3 * (2 * a - b);
+        auto eqnE = (2 * a - b) + (2 * a - b);
+        auto eqnF = (2 * a - b) - (2 * a - b);
+
+        REQUIRE(eqnA.size() == 2);
+        REQUIRE(eqnB.size() == 2);
+        REQUIRE(eqnC.size() == 2);
+
+        REQUIRE(getField(eqnA.explicitOperation(size)) == 4);
+        REQUIRE(getField(eqnB.explicitOperation(size)) == 12);
+        REQUIRE(getField(eqnC.explicitOperation(size)) == 2);
+        REQUIRE(getField(eqnD.explicitOperation(size)) == 6);
+        REQUIRE(getField(eqnE.explicitOperation(size)) == 4);
+        REQUIRE(getField(eqnF.explicitOperation(size)) == 0);
+    }
+}
diff --git a/test/dsl/operator.cpp b/test/dsl/operator.cpp
index b48c4a1ee..5203aa9c9 100644
--- a/test/dsl/operator.cpp
+++ b/test/dsl/operator.cpp
@@ -23,7 +23,7 @@ TEST_CASE("Operator")
 
         std::vector<fvcc::VolumeBoundary<NeoFOAM::scalar>> bcs {};
         auto vf = VolumeField(exec, mesh, fA, bf, bcs);
-        auto b = Dummy(exec, vf);
+        auto b = Dummy(vf);
 
         REQUIRE(b.getName() == "Dummy");
         REQUIRE(b.getType() == Operator::Type::Explicit);
@@ -38,13 +38,13 @@ TEST_CASE("Operator")
         BoundaryFields bf(exec, mesh.nBoundaryFaces(), mesh.nBoundaries());
         auto vf = VolumeField(exec, mesh, fA, bf, bcs);
 
-        auto c = 2 * Dummy(exec, vf);
-        auto d = fB * Dummy(exec, vf);
-        auto e = Coeff(-3, fB) * Dummy(exec, vf);
+        auto c = 2 * Dummy(vf);
+        auto d = fB * Dummy(vf);
+        auto e = Coeff(-3, fB) * Dummy(vf);
 
-        auto coeffc = c.getCoefficient();
-        auto coeffd = d.getCoefficient();
-        auto coeffE = e.getCoefficient();
+        [[maybe_unused]] auto coeffC = c.getCoefficient();
+        [[maybe_unused]] auto coeffD = d.getCoefficient();
+        [[maybe_unused]] auto coeffE = e.getCoefficient();
 
         Field source(exec, 1, 2.0);
         c.explicitOperation(source);
diff --git a/test/finiteVolume/cellCentred/timeIntegration/CMakeLists.txt b/test/dsl/timeIntegration/CMakeLists.txt
similarity index 100%
rename from test/finiteVolume/cellCentred/timeIntegration/CMakeLists.txt
rename to test/dsl/timeIntegration/CMakeLists.txt
diff --git a/test/dsl/timeIntegration/timeIntegration.cpp b/test/dsl/timeIntegration/timeIntegration.cpp
new file mode 100644
index 000000000..4abb44e66
--- /dev/null
+++ b/test/dsl/timeIntegration/timeIntegration.cpp
@@ -0,0 +1,55 @@
+// SPDX-License-Identifier: MIT
+// SPDX-FileCopyrightText: 2024 NeoFOAM authors
+
+#define CATCH_CONFIG_RUNNER // Define this before including catch.hpp to create
+                            // a custom main
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/generators/catch_generators_all.hpp>
+
+#include "../common.hpp"
+
+#include "NeoFOAM/core/dictionary.hpp"
+#include "NeoFOAM/core/parallelAlgorithms.hpp"
+#include "NeoFOAM/timeIntegration/forwardEuler.hpp"
+#include "NeoFOAM/dsl/solver.hpp"
+#include "NeoFOAM/dsl/ddt.hpp"
+
+TEST_CASE("TimeIntegration")
+{
+    NeoFOAM::Executor exec = GENERATE(
+        NeoFOAM::Executor(NeoFOAM::SerialExecutor {}),
+        NeoFOAM::Executor(NeoFOAM::CPUExecutor {}),
+        NeoFOAM::Executor(NeoFOAM::GPUExecutor {})
+    );
+
+    std::string execName = std::visit([](auto e) { return e.print(); }, exec);
+    auto mesh = NeoFOAM::createSingleCellMesh(exec);
+
+    NeoFOAM::Dictionary fvSchemes;
+    NeoFOAM::Dictionary ddtSchemes;
+    ddtSchemes.insert("type", std::string("forwardEuler"));
+    fvSchemes.insert("ddtSchemes", ddtSchemes);
+    NeoFOAM::Dictionary fvSolution;
+
+    Field fA(exec, 1, 2.0);
+    BoundaryFields bf(exec, mesh.nBoundaryFaces(), mesh.nBoundaries());
+
+    std::vector<fvcc::VolumeBoundary<NeoFOAM::scalar>> bcs {};
+    auto vf = VolumeField(exec, mesh, fA, bf, bcs);
+
+    SECTION("Create expression and perform explicitOperation on " + execName)
+    {
+        auto dummy = Dummy(vf);
+        Operator ddtOperator = NeoFOAM::dsl::temporal::ddt(vf);
+
+        // ddt(U) = f
+        auto eqn = ddtOperator + dummy;
+        double dt {2.0};
+
+        // int(ddt(U)) + f = 0
+        // (U^1-U^0)/dt = -f
+        // U^1 = - f * dt + U^0, where dt = 2, f=1, U^0=2.0 -> U^1=-2.0
+        solve(eqn, vf, dt, fvSchemes, fvSolution);
+        REQUIRE(getField(vf.internalField()) == -2.0);
+    }
+}
diff --git a/test/finiteVolume/CMakeLists.txt b/test/finiteVolume/CMakeLists.txt
index 7213b6f5a..564d38f7b 100644
--- a/test/finiteVolume/CMakeLists.txt
+++ b/test/finiteVolume/CMakeLists.txt
@@ -4,4 +4,3 @@
 add_subdirectory(cellCentred/fields)
 add_subdirectory(cellCentred/boundary)
 add_subdirectory(cellCentred/interpolation)
-add_subdirectory(cellCentred/timeIntegration)
diff --git a/test/finiteVolume/cellCentred/timeIntegration/timeIntegration.cpp b/test/finiteVolume/cellCentred/timeIntegration/timeIntegration.cpp
deleted file mode 100644
index 95a259cae..000000000
--- a/test/finiteVolume/cellCentred/timeIntegration/timeIntegration.cpp
+++ /dev/null
@@ -1,94 +0,0 @@
-// SPDX-License-Identifier: MIT
-// SPDX-FileCopyrightText: 2024 NeoFOAM authors
-
-#define CATCH_CONFIG_RUNNER // Define this before including catch.hpp to create
-                            // a custom main
-#include <catch2/catch_session.hpp>
-#include <catch2/catch_test_macros.hpp>
-#include <catch2/generators/catch_generators_all.hpp>
-
-#include "NeoFOAM/core/dictionary.hpp"
-#include "NeoFOAM/core/parallelAlgorithms.hpp"
-#include "NeoFOAM/finiteVolume/cellCentred/timeIntegration/timeIntegration.hpp"
-
-namespace dsl = NeoFOAM::DSL;
-
-class Divergence : public NeoFOAM::DSL::OperatorMixin
-{
-
-public:
-
-    std::string display() const { return "Divergence"; }
-
-    void explicitOperation(NeoFOAM::Field<NeoFOAM::scalar>& source, NeoFOAM::scalar scale)
-    {
-        auto sourceField = source.span();
-        NeoFOAM::parallelFor(
-            source.exec(),
-            {0, source.size()},
-            KOKKOS_LAMBDA(const size_t i) { sourceField[i] += 1.0 * scale; }
-        );
-    }
-
-    dsl::Operator::Type getType() const { return termType_; }
-
-    fvcc::VolumeField<NeoFOAM::scalar>* volumeField() const { return nullptr; }
-
-    const NeoFOAM::Executor& exec() const { return exec_; }
-
-    std::size_t nCells() const { return nCells_; }
-
-    dsl::Operator::Type termType_;
-
-    NeoFOAM::Executor exec_;
-
-    std::size_t nCells_;
-};
-
-class TimeOperator
-{
-
-public:
-
-    std::string display() const { return "TimeOperator"; }
-
-    void explicitOperation(NeoFOAM::Field<NeoFOAM::scalar>& source, NeoFOAM::scalar scale)
-    {
-        auto sourceField = source.span();
-        NeoFOAM::parallelFor(
-            source.exec(),
-            {0, source.size()},
-            KOKKOS_LAMBDA(const size_t i) { sourceField[i] += 1 * scale; }
-        );
-    }
-
-    dsl::Operator::Type getType() const { return termType_; }
-
-    fvcc::VolumeField<NeoFOAM::scalar>* volumeField() const { return nullptr; }
-
-    const NeoFOAM::Executor& exec() const { return exec_; }
-
-    std::size_t nCells() const { return nCells_; }
-
-    dsl::Operator::Type termType_;
-    const NeoFOAM::Executor exec_;
-    std::size_t nCells_;
-};
-
-
-TEST_CASE("TimeIntegration")
-{
-    auto exec = NeoFOAM::SerialExecutor();
-    namespace fvcc = NeoFOAM::finiteVolume::cellCentred;
-
-    NeoFOAM::Dictionary dict;
-    dict.insert("type", std::string("forwardEuler"));
-
-    // dsl::Operator divOperator = Divergence(dsl::Operator::Type::Explicit, exec, 1);
-
-    // dsl::Operator ddtOperator = TimeOperator(dsl::Operator::Type::Temporal, exec, 1);
-
-    // dsl::EqnSystem eqnSys = ddtOperator + divOperator;
-
-    // fvcc::TimeIntegration timeIntergrator(eqnSys, dict);
-}