diff --git a/CHANGELOG.md b/CHANGELOG.md
index 97a54c1ab..d92d41391 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,6 @@
 # Version 0.2.0 (unreleased)
 ## Features
+- support for implicit operators in the DSL [#246](https://github.com/exasim-project/NeoFOAM/pull/246)
 - add segmentedField to represent vector of vectors [#202](https://github.com/exasim-project/NeoFOAM/pull/202)
 - Adds a minimal implementation linear algebra functionality [#219](https://github.com/exasim-project/NeoFOAM/pull/219)
 ## Fixes
diff --git a/include/NeoFOAM/dsl/expression.hpp b/include/NeoFOAM/dsl/expression.hpp
index da84754d5..749964968 100644
--- a/include/NeoFOAM/dsl/expression.hpp
+++ b/include/NeoFOAM/dsl/expression.hpp
@@ -8,9 +8,12 @@
 
 #include "NeoFOAM/core/primitives/scalar.hpp"
 #include "NeoFOAM/fields/field.hpp"
+#include "NeoFOAM/linearAlgebra.hpp"
 #include "NeoFOAM/dsl/operator.hpp"
 #include "NeoFOAM/core/error.hpp"
 
+namespace la = NeoFOAM::la;
+
 namespace NeoFOAM::dsl
 {
 
@@ -61,6 +64,21 @@ class Expression
         return source;
     }
 
+    /* @brief perform all implicit operation and accumulate the result */
+    la::LinearSystem<scalar, localIdx> implicitOperation()
+    {
+        if (implicitOperators_.empty())
+        {
+            NF_ERROR_EXIT("No implicit operators in the expression");
+        }
+        auto ls = implicitOperators_[0].createEmptyLinearSystem();
+        for (auto& oper : implicitOperators_)
+        {
+            oper.implicitOperation(ls);
+        }
+        return ls;
+    }
+
     void addOperator(const Operator& oper)
     {
         switch (oper.getType())
diff --git a/include/NeoFOAM/dsl/operator.hpp b/include/NeoFOAM/dsl/operator.hpp
index ef3e12fb8..a5693ee30 100644
--- a/include/NeoFOAM/dsl/operator.hpp
+++ b/include/NeoFOAM/dsl/operator.hpp
@@ -7,9 +7,12 @@
 
 #include "NeoFOAM/core/primitives/scalar.hpp"
 #include "NeoFOAM/fields/field.hpp"
+#include "NeoFOAM/linearAlgebra.hpp"
 #include "NeoFOAM/core/input.hpp"
 #include "NeoFOAM/dsl/coeff.hpp"
 
+namespace la = NeoFOAM::la;
+
 namespace NeoFOAM::dsl
 {
 
@@ -27,6 +30,13 @@ concept HasExplicitOperator = requires(T t) {
     } -> std::same_as<void>; // Adjust return type and arguments as needed
 };
 
+template<typename T>
+concept HasImplicitOperator = requires(T t) {
+    {
+        t.implicitOperation(std::declval<la::LinearSystem<scalar, localIdx>&>())
+    } -> std::same_as<void>; // Adjust return type and arguments as needed
+};
+
 /* @class Operator
  * @brief A class to represent an operator in NeoFOAMs dsl
  *
@@ -62,6 +72,10 @@ class Operator
 
     void temporalOperation(Field<scalar>& field);
 
+    void implicitOperation(la::LinearSystem<scalar, localIdx>& ls);
+
+    la::LinearSystem<scalar, localIdx> createEmptyLinearSystem() const;
+
     /* returns the fundamental type of an operator, ie explicit, implicit, temporal */
     Operator::Type getType() const;
 
@@ -89,6 +103,10 @@ class Operator
 
         virtual void temporalOperation(Field<scalar>& field) = 0;
 
+        virtual void implicitOperation(la::LinearSystem<scalar, localIdx>& ls) = 0;
+
+        virtual la::LinearSystem<scalar, localIdx> createEmptyLinearSystem() const = 0;
+
         /* @brief Given an input this function reads required coeffs */
         virtual void build(const Input& input) = 0;
 
@@ -137,6 +155,35 @@ class Operator
             }
         }
 
+        virtual void implicitOperation(la::LinearSystem<scalar, localIdx>& ls) override
+        {
+            if constexpr (HasImplicitOperator<ConcreteOperatorType>)
+            {
+                concreteOp_.implicitOperation(ls);
+            }
+        }
+
+        virtual la::LinearSystem<scalar, localIdx> createEmptyLinearSystem() const override
+        {
+            if constexpr (HasImplicitOperator<ConcreteOperatorType>)
+            {
+                return concreteOp_.createEmptyLinearSystem();
+            }
+            throw std::runtime_error("Implicit operation not implemented");
+            // only need to avoid compiler warning about missing return statement
+            // this code path should never be reached as we call implicitOperation on an explicit
+            // operator
+            NeoFOAM::Field<NeoFOAM::scalar> values(exec(), 1, 0.0);
+            NeoFOAM::Field<NeoFOAM::localIdx> colIdx(exec(), 1, 0);
+            NeoFOAM::Field<NeoFOAM::localIdx> rowPtrs(exec(), 2, 0);
+            NeoFOAM::la::CSRMatrix<NeoFOAM::scalar, NeoFOAM::localIdx> csrMatrix(
+                values, colIdx, rowPtrs
+            );
+
+            NeoFOAM::Field<NeoFOAM::scalar> rhs(exec(), 1, 0.0);
+            return la::LinearSystem<scalar, localIdx>(csrMatrix, rhs, "custom");
+        }
+
         /* @brief Given an input this function reads required coeffs */
         virtual void build(const Input& input) override { concreteOp_.build(input); }
 
diff --git a/include/NeoFOAM/linearAlgebra/CSRMatrix.hpp b/include/NeoFOAM/linearAlgebra/CSRMatrix.hpp
index 9ff89e527..bbdb9204f 100644
--- a/include/NeoFOAM/linearAlgebra/CSRMatrix.hpp
+++ b/include/NeoFOAM/linearAlgebra/CSRMatrix.hpp
@@ -37,6 +37,11 @@ class CSRMatrix
 
     [[nodiscard]] std::size_t nColIdxs() const { return colIdxs_.size(); }
 
+    [[nodiscard]] CSRMatrix copyToHost() const
+    {
+        return CSRMatrix(values_.copyToHost(), colIdxs_.copyToHost(), rowPtrs_.copyToHost());
+    }
+
     [[nodiscard]] std::span<ValueType> values() { return values_.span(); }
     [[nodiscard]] std::span<IndexType> colIdxs() { return colIdxs_.span(); }
     [[nodiscard]] std::span<IndexType> rowPtrs() { return rowPtrs_.span(); }
diff --git a/include/NeoFOAM/linearAlgebra/linearSystem.hpp b/include/NeoFOAM/linearAlgebra/linearSystem.hpp
index 5be2654aa..38d16799a 100644
--- a/include/NeoFOAM/linearAlgebra/linearSystem.hpp
+++ b/include/NeoFOAM/linearAlgebra/linearSystem.hpp
@@ -22,8 +22,12 @@ class LinearSystem
 {
 public:
 
-    LinearSystem(const CSRMatrix<ValueType, IndexType>& matrix, const Field<ValueType>& rhs)
-        : matrix_(matrix), rhs_(rhs)
+    LinearSystem(
+        const CSRMatrix<ValueType, IndexType>& matrix,
+        const Field<ValueType>& rhs,
+        const std::string& sparsityPattern
+    )
+        : matrix_(matrix), rhs_(rhs), sparsityPattern_(sparsityPattern)
     {
         NF_ASSERT(matrix.exec() == rhs.exec(), "Executors are not the same");
         NF_ASSERT(matrix.nRows() == rhs.size(), "Matrix and RHS size mismatch");
@@ -41,12 +45,20 @@ class LinearSystem
     [[nodiscard]] const CSRMatrix<ValueType, IndexType>& matrix() const { return matrix_; }
     [[nodiscard]] const Field<ValueType>& rhs() const { return rhs_; }
 
+    [[nodiscard]] std::string sparsityPattern() const { return sparsityPattern_; }
+
+    [[nodiscard]] LinearSystem copyToHost() const
+    {
+        return LinearSystem(matrix_.copyToHost(), rhs_.copyToHost(), sparsityPattern_);
+    }
+
     const Executor& exec() const { return matrix_.exec(); }
 
 private:
 
     CSRMatrix<ValueType, IndexType> matrix_;
     Field<ValueType> rhs_;
+    std::string sparsityPattern_;
 };
 
 
diff --git a/src/dsl/operator.cpp b/src/dsl/operator.cpp
index b6be5ace1..b85c018d9 100644
--- a/src/dsl/operator.cpp
+++ b/src/dsl/operator.cpp
@@ -15,6 +15,16 @@ void Operator::explicitOperation(Field<scalar>& source) { model_->explicitOperat
 
 void Operator::temporalOperation(Field<scalar>& field) { model_->temporalOperation(field); }
 
+void Operator::implicitOperation(la::LinearSystem<scalar, localIdx>& ls)
+{
+    model_->implicitOperation(ls);
+}
+
+la::LinearSystem<scalar, localIdx> Operator::createEmptyLinearSystem() const
+{
+    return model_->createEmptyLinearSystem();
+}
+
 Operator::Type Operator::getType() const { return model_->getType(); }
 
 Coeff& Operator::getCoefficient() { return model_->getCoefficient(); }
diff --git a/test/dsl/common.hpp b/test/dsl/common.hpp
index c2380af35..7fc3e0ba1 100644
--- a/test/dsl/common.hpp
+++ b/test/dsl/common.hpp
@@ -28,6 +28,8 @@ class Dummy : public OperatorMixin
 
     Dummy(VolumeField& field) : OperatorMixin(field.exec(), field, Operator::Type::Explicit) {}
 
+    Dummy(VolumeField& field, Operator::Type type) : OperatorMixin(field.exec(), field, type) {}
+
     void explicitOperation(Field& source)
     {
         auto sourceSpan = source.span();
@@ -40,6 +42,44 @@ class Dummy : public OperatorMixin
         );
     }
 
+    void implicitOperation(la::LinearSystem<NeoFOAM::scalar, NeoFOAM::localIdx>& ls)
+    {
+        auto values = ls.matrix().values();
+        auto rhs = ls.rhs().span();
+        auto fieldSpan = field_.internalField().span();
+        auto coeff = getCoefficient();
+
+        // update diag
+        NeoFOAM::parallelFor(
+            exec(),
+            {0, values.size()},
+            KOKKOS_LAMBDA(const size_t i) { values[i] += coeff[i] * fieldSpan[i]; }
+        );
+
+        // update rhs
+        NeoFOAM::parallelFor(
+            exec(),
+            ls.rhs().range(),
+            KOKKOS_LAMBDA(const size_t i) { rhs[i] += coeff[i] * fieldSpan[i]; }
+        );
+    }
+
+    la::LinearSystem<NeoFOAM::scalar, NeoFOAM::localIdx> createEmptyLinearSystem() const
+    {
+        NeoFOAM::Field<NeoFOAM::scalar> values(exec(), 1, 0.0);
+        NeoFOAM::Field<NeoFOAM::localIdx> colIdx(exec(), 1, 0.0);
+        NeoFOAM::Field<NeoFOAM::localIdx> rowPtrs(exec(), {0, 1});
+        NeoFOAM::la::CSRMatrix<NeoFOAM::scalar, NeoFOAM::localIdx> csrMatrix(
+            values, colIdx, rowPtrs
+        );
+
+        NeoFOAM::Field<NeoFOAM::scalar> rhs(exec(), 1, 0.0);
+        NeoFOAM::la::LinearSystem<NeoFOAM::scalar, NeoFOAM::localIdx> linearSystem(
+            csrMatrix, rhs, "diagonal"
+        );
+        return linearSystem;
+    }
+
     std::string getName() const { return "Dummy"; }
 };
 
@@ -48,3 +88,16 @@ NeoFOAM::scalar getField(const NeoFOAM::Field<NeoFOAM::scalar>& source)
     auto sourceField = source.copyToHost();
     return sourceField.span()[0];
 }
+
+
+NeoFOAM::scalar getDiag(const la::LinearSystem<NeoFOAM::scalar, NeoFOAM::localIdx> ls)
+{
+    auto hostLs = ls.copyToHost();
+    return hostLs.matrix().values()[0];
+}
+
+NeoFOAM::scalar getRhs(const la::LinearSystem<NeoFOAM::scalar, NeoFOAM::localIdx> ls)
+{
+    auto hostLs = ls.copyToHost();
+    return hostLs.rhs().span()[0];
+}
diff --git a/test/dsl/expression.cpp b/test/dsl/expression.cpp
index 4d5a31eea..d7a4a913f 100644
--- a/test/dsl/expression.cpp
+++ b/test/dsl/expression.cpp
@@ -26,11 +26,12 @@ TEST_CASE("Expression")
     auto vf = VolumeField(exec, "vf", 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 a = Dummy(vf);
+        auto b = Dummy(vf);
+
         auto eqnA = a + b;
         auto eqnB = fB * Dummy(vf) + 2 * Dummy(vf);
         auto eqnC = Expression(2 * a - b);
@@ -49,4 +50,40 @@ TEST_CASE("Expression")
         REQUIRE(getField(eqnE.explicitOperation(size)) == 4);
         REQUIRE(getField(eqnF.explicitOperation(size)) == 0);
     }
+
+    SECTION("Create equation and perform implicitOperation on " + execName)
+    {
+        auto a = Dummy(vf, Operator::Type::Implicit);
+        auto b = Dummy(vf, Operator::Type::Implicit);
+
+        auto eqnA = a + b;
+        auto eqnB =
+            fB * Dummy(vf, Operator::Type::Implicit) + 2 * Dummy(vf, Operator::Type::Implicit);
+        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(getDiag(eqnA.implicitOperation()) == 4);
+        REQUIRE(getRhs(eqnA.implicitOperation()) == 4);
+
+        REQUIRE(getDiag(eqnB.implicitOperation()) == 12);
+        REQUIRE(getRhs(eqnB.implicitOperation()) == 12);
+
+        REQUIRE(getDiag(eqnC.implicitOperation()) == 2);
+        REQUIRE(getRhs(eqnC.implicitOperation()) == 2);
+
+        REQUIRE(getDiag(eqnD.implicitOperation()) == 6);
+        REQUIRE(getRhs(eqnD.implicitOperation()) == 6);
+
+        REQUIRE(getDiag(eqnE.implicitOperation()) == 4);
+        REQUIRE(getRhs(eqnE.implicitOperation()) == 4);
+
+        REQUIRE(getDiag(eqnF.implicitOperation()) == 0);
+        REQUIRE(getRhs(eqnF.implicitOperation()) == 0);
+    }
 }
diff --git a/test/dsl/operator.cpp b/test/dsl/operator.cpp
index 304fa9ff8..1d2f1ed31 100644
--- a/test/dsl/operator.cpp
+++ b/test/dsl/operator.cpp
@@ -29,7 +29,7 @@ TEST_CASE("Operator")
         REQUIRE(b.getType() == Operator::Type::Explicit);
     }
 
-    SECTION("Supports Coefficients" + execName)
+    SECTION("Supports Coefficients Explicit " + execName)
     {
         std::vector<fvcc::VolumeBoundary<NeoFOAM::scalar>> bcs {};
 
@@ -63,4 +63,68 @@ TEST_CASE("Operator")
         auto hostSourceE = source.copyToHost();
         REQUIRE(hostSourceE.span()[0] == -2.0);
     }
+
+    SECTION("Implicit Operations " + execName)
+    {
+        Field fA(exec, 1, 2.0);
+        BoundaryFields bf(exec, mesh.nBoundaryFaces(), mesh.nBoundaries());
+
+        std::vector<fvcc::VolumeBoundary<NeoFOAM::scalar>> bcs {};
+        auto vf = VolumeField(exec, "vf", mesh, fA, bf, bcs);
+        auto b = Dummy(vf, Operator::Type::Implicit);
+
+        REQUIRE(b.getName() == "Dummy");
+        REQUIRE(b.getType() == Operator::Type::Implicit);
+
+        auto ls = b.createEmptyLinearSystem();
+        REQUIRE(ls.matrix().nValues() == 1);
+        REQUIRE(ls.matrix().nColIdxs() == 1);
+        REQUIRE(ls.matrix().nRows() == 1);
+    }
+
+    SECTION("Supports Coefficients Implicit " + execName)
+    {
+        std::vector<fvcc::VolumeBoundary<NeoFOAM::scalar>> bcs {};
+
+        Field fA(exec, 1, 2.0);
+        Field fB(exec, 1, 2.0);
+        BoundaryFields bf(exec, mesh.nBoundaryFaces(), mesh.nBoundaries());
+        auto vf = VolumeField(exec, "vf", mesh, fA, bf, bcs);
+
+        auto c = 2 * Dummy(vf, Operator::Type::Implicit);
+        auto d = fB * Dummy(vf, Operator::Type::Implicit);
+        auto e = Coeff(-3, fB) * Dummy(vf, Operator::Type::Implicit);
+
+        [[maybe_unused]] auto coeffC = c.getCoefficient();
+        [[maybe_unused]] auto coeffD = d.getCoefficient();
+        [[maybe_unused]] auto coeffE = e.getCoefficient();
+
+        // Field source(exec, 1, 2.0);
+        auto ls = c.createEmptyLinearSystem();
+        c.implicitOperation(ls);
+
+        // c = 2 * 2
+        auto hostRhsC = ls.rhs().copyToHost();
+        REQUIRE(hostRhsC.span()[0] == 4.0);
+        auto hostLsC = ls.copyToHost();
+        REQUIRE(hostLsC.matrix().values()[0] == 4.0);
+
+
+        // d= 2 * 2
+        ls = d.createEmptyLinearSystem();
+        d.implicitOperation(ls);
+        auto hostRhsD = ls.rhs().copyToHost();
+        REQUIRE(hostRhsD.span()[0] == 4.0);
+        auto hostLsD = ls.copyToHost();
+        REQUIRE(hostLsD.matrix().values()[0] == 4.0);
+
+
+        // e = - -3 * 2 * 2 = -12
+        ls = e.createEmptyLinearSystem();
+        e.implicitOperation(ls);
+        auto hostRhsE = ls.rhs().copyToHost();
+        REQUIRE(hostRhsE.span()[0] == -12.0);
+        auto hostLsE = ls.copyToHost();
+        REQUIRE(hostLsE.matrix().values()[0] == -12.0);
+    }
 }
diff --git a/test/linearAlgebra/linearSystem.cpp b/test/linearAlgebra/linearSystem.cpp
index 0b491e9c9..f3ecfe977 100644
--- a/test/linearAlgebra/linearSystem.cpp
+++ b/test/linearAlgebra/linearSystem.cpp
@@ -31,7 +31,9 @@ TEST_CASE("LinearSystem")
         );
 
         NeoFOAM::Field<NeoFOAM::scalar> rhs(exec, 3, 0.0);
-        NeoFOAM::la::LinearSystem<NeoFOAM::scalar, NeoFOAM::localIdx> linearSystem(csrMatrix, rhs);
+        NeoFOAM::la::LinearSystem<NeoFOAM::scalar, NeoFOAM::localIdx> linearSystem(
+            csrMatrix, rhs, "custom"
+        );
 
         REQUIRE(linearSystem.matrix().values().size() == 9);
         REQUIRE(linearSystem.matrix().nValues() == 9);
@@ -54,7 +56,9 @@ TEST_CASE("LinearSystem")
         );
 
         NeoFOAM::Field<NeoFOAM::scalar> rhs(exec, 3, 0.0);
-        NeoFOAM::la::LinearSystem<NeoFOAM::scalar, NeoFOAM::localIdx> linearSystem(csrMatrix, rhs);
+        NeoFOAM::la::LinearSystem<NeoFOAM::scalar, NeoFOAM::localIdx> linearSystem(
+            csrMatrix, rhs, "testing"
+        );
         NeoFOAM::Field<NeoFOAM::scalar> x(exec, {1.0, 2.0, 3.0});
 
         NeoFOAM::Field<NeoFOAM::scalar> y = NeoFOAM::la::SpMV(linearSystem, x);
@@ -68,7 +72,7 @@ TEST_CASE("LinearSystem")
         // test with non-zero rhs
         NeoFOAM::Field<NeoFOAM::scalar> rhs2(exec, {1.0, 2.0, 3.0});
         NeoFOAM::la::LinearSystem<NeoFOAM::scalar, NeoFOAM::localIdx> linearSystem2(
-            csrMatrix, rhs2
+            csrMatrix, rhs2, "testing"
         );
         y = NeoFOAM::la::SpMV(linearSystem2, x);
         yHost = y.copyToHost();