diff --git a/framework/doc/content/source/mfem/equation_systems/EquationSystem.md b/framework/doc/content/source/mfem/equation_systems/EquationSystem.md
index 5e3858cf3075..55a8819df153 100644
--- a/framework/doc/content/source/mfem/equation_systems/EquationSystem.md
+++ b/framework/doc/content/source/mfem/equation_systems/EquationSystem.md
@@ -32,6 +32,18 @@ form. [`mfem::BilinearFormIntegrators`](https://mfem.org/bilininteg/) add contri
 $A_{ij}(\varphi_i, \phi_j)$ and [`mfem::LinearFormIntegrators`](https://mfem.org/lininteg/) add
 contributions to $b_i(\varphi_i)$ when assembled.
 
+[MFEMKernels](source/mfem/kernels/MFEMKernel.md) can also contribute to domain integrators for 
+non-linear actions. This allows to form the residual $\mathcal{L}(u)$ for non-linear Newton's 
+methood as shown below
+
+!equation
+{\mathbf{J}\left(\vec{u}_n\right) \delta \vec{u}_{n+1}=-\vec{R}\left(\vec{u}_n\right)}
+
+!equation
+{\vec{u}_{n+1}=\vec{u}_n+\delta \vec{u}_{n+1}}
+
+where $\mathbf{J}$ is the Jacobian, and $\delta \vec{u}$ is the incremental solution.
+
 !if-end!
 
 !else
diff --git a/framework/doc/content/source/mfem/kernels/MFEMKernel.md b/framework/doc/content/source/mfem/kernels/MFEMKernel.md
index 3eabf1e5fd77..9305f09c0884 100644
--- a/framework/doc/content/source/mfem/kernels/MFEMKernel.md
+++ b/framework/doc/content/source/mfem/kernels/MFEMKernel.md
@@ -22,7 +22,9 @@ variable names is the same as the set of trial variable names for a square syste
 
 `MFEMKernel` is a purely virtual base class. Derived classes should override the `createBFIntegrator`
 and/or the `createLFIntegrator` methods to return a `BilinearFormIntegrator` and/or a
-`LinearFormIntegrator` (respectively) to add to the `EquationSystem`.
+`LinearFormIntegrator` (respectively) to add to the `EquationSystem`. Derived classes could also 
+override `createNLAIntegrator` to solve non-linear problems using 
+[EquationSystem](source/mfem/equation_systems/EquationSystem.md).
 
 !if-end!
 
diff --git a/framework/doc/content/source/mfem/kernels/MFEMNLDiffusionKernel.md b/framework/doc/content/source/mfem/kernels/MFEMNLDiffusionKernel.md
new file mode 100644
index 000000000000..35433988528e
--- /dev/null
+++ b/framework/doc/content/source/mfem/kernels/MFEMNLDiffusionKernel.md
@@ -0,0 +1,37 @@
+# MFEMNLDiffusionKernel
+
+!if! function=hasCapability('mfem')
+
+## Overview
+
+Adds the domain integrator for integrating the non-linear action
+
+!equation
+(k(u)\vec\nabla v, \vec\nabla v)_\Omega \,\,\, \forall v \in V
+
+Adds the domain integrator for integrating the bilinear form
+
+!equation
+(k(u)\vec\nabla v, \vec\nabla v)_\Omega + (k'(u) v, \vec\nabla u \vec\nabla v)_\Omega \,\,\, \forall v \in V
+
+where $u, v \in H^1$ and $k(u)$ is a scalar non-linear diffusivity coefficient.
+
+The above terms arises from the weak form of the non-linear operator
+
+!equation
+- \vec\nabla \cdot \left( k(u) \vec\nabla u \right)
+
+## Example Input File Syntax
+
+!listing mfem/kernels/nldiffusion.i block=/Kernels
+
+!syntax parameters /Kernels/MFEMNLDiffusionKernel
+
+!syntax inputs /Kernels/MFEMNLDiffusionKernel
+
+!syntax children /Kernels/MFEMNLDiffusionKernel
+
+!if-end!
+
+!else
+!include mfem/mfem_warning.md
diff --git a/framework/include/mfem/bcs/MFEMIntegratedBC.h b/framework/include/mfem/bcs/MFEMIntegratedBC.h
index 954468fbfa4a..4358c082d98e 100644
--- a/framework/include/mfem/bcs/MFEMIntegratedBC.h
+++ b/framework/include/mfem/bcs/MFEMIntegratedBC.h
@@ -24,6 +24,9 @@ class MFEMIntegratedBC : public MFEMBoundaryCondition
   /// Create MFEM integrator to apply to the RHS of the weak form. Ownership managed by the caller.
   virtual mfem::LinearFormIntegrator * createLFIntegrator() = 0;
 
+  /// Create MFEM integrator to apply non-linear residual form. Ownership managed by the caller.
+  virtual mfem::LinearFormIntegrator * createNLAIntegrator() { return nullptr; };
+
   /// Create MFEM integrator to apply to the LHS of the weak form. Ownership managed by the caller.
   virtual mfem::BilinearFormIntegrator * createBFIntegrator() = 0;
 
diff --git a/framework/include/mfem/equation_systems/ComplexEquationSystem.h b/framework/include/mfem/equation_systems/ComplexEquationSystem.h
index 01688b60320d..a3f810000e7c 100644
--- a/framework/include/mfem/equation_systems/ComplexEquationSystem.h
+++ b/framework/include/mfem/equation_systems/ComplexEquationSystem.h
@@ -28,6 +28,9 @@ class ComplexEquationSystem : public EquationSystem
                     ComplexGridFunctions & cmplx_gridfunctions,
                     mfem::AssemblyLevel assembly_level) override;
 
+  ///Nonlinear Mult (Used by Newton-solver not necessarily nonlinear)
+  virtual void Mult(const mfem::Vector & x, mfem::Vector & y) const override;
+
   /// Build all forms comprising this EquationSystem
   virtual void BuildEquationSystem() override;
 
diff --git a/framework/include/mfem/equation_systems/EquationSystem.h b/framework/include/mfem/equation_systems/EquationSystem.h
index 531a18d6033f..2f669e06a34a 100644
--- a/framework/include/mfem/equation_systems/EquationSystem.h
+++ b/framework/include/mfem/equation_systems/EquationSystem.h
@@ -5,7 +5,7 @@
 //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
 //*
 //* Licensed under LGPL 2.1, please see LICENSE for details
-//* https://www.gnu.org/licenses/lgpl-2.1.html
+//* https://www.gnu.org/licenses/lgpl-2.1.html*/
 
 #ifdef MOOSE_MFEM_ENABLED
 
@@ -85,17 +85,21 @@ class EquationSystem : public mfem::Operator
                                 mfem::Array<int> & global_ess_markers);
   /// Update all essentially constrained true DoF markers and values on boundaries
   virtual void ApplyEssentialBCs();
-
   /// Perform trivial eliminations of coupled variables lacking corresponding test variables
   virtual void EliminateCoupledVariables();
   /// Build linear forms and eliminate constrained DoFs
   virtual void BuildLinearForms();
+  /// Build non-linear action forms
+  virtual void BuildNonLinearActionForms();
   /// Build bilinear forms (diagonal Jacobian contributions)
   virtual void BuildBilinearForms();
   /// Build mixed bilinear forms (off-diagonal Jacobian contributions)
   virtual void BuildMixedBilinearForms();
   /// Build all forms comprising this EquationSystem
   virtual void BuildEquationSystem();
+  // Reassemble the Jacobian matrix/Operator (Assuming Gridfunction dependant Operator)
+  virtual void ReassembleJacobian(mfem::BlockVector & x, mfem::BlockVector & rhs);
+  virtual void UpdateJacobian() const;
 
   /// Form linear system and jacobian operator based on on- and off-diagonal bilinear form
   /// contributions, populate solution and RHS vectors of true DoFs, and apply constraints.
@@ -130,6 +134,11 @@ class EquationSystem : public mfem::Operator
       std::shared_ptr<mfem::ParLinearForm> form,
       NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map);
 
+  void ApplyDomainNLAFIntegrators(
+      const std::string & test_var_name,
+      std::shared_ptr<mfem::ParLinearForm> form,
+      NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map);
+
   template <class FormType>
   void ApplyBoundaryBLFIntegrators(
       const std::string & trial_var_name,
@@ -145,6 +154,12 @@ class EquationSystem : public mfem::Operator
       NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
           integrated_bc_map);
 
+  void ApplyBoundaryNLAFIntegrators(
+      const std::string & test_var_name,
+      std::shared_ptr<mfem::ParLinearForm> form,
+      NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
+          integrated_bc_map);
+
   /// Names of all trial variables of kernels and boundary conditions
   /// added to this EquationSystem.
   std::vector<std::string> _coupled_var_names;
@@ -167,8 +182,8 @@ class EquationSystem : public mfem::Operator
   NamedFieldsMap<mfem::ParBilinearForm> _blfs;
   NamedFieldsMap<mfem::ParLinearForm> _lfs;
   NamedFieldsMap<mfem::ParNonlinearForm> _nlfs;
-  NamedFieldsMap<NamedFieldsMap<mfem::ParMixedBilinearForm>>
-      _mblfs; // named according to trial variable
+  NamedFieldsMap<mfem::ParLinearForm> _nlAs;
+  NamedFieldsMap<NamedFieldsMap<mfem::ParMixedBilinearForm>> _mblfs; // named according to trial var
 
   /// Gridfunctions holding essential constraints from Dirichlet BCs
   std::vector<std::unique_ptr<mfem::ParGridFunction>> _var_ess_constraints;
@@ -185,10 +200,18 @@ class EquationSystem : public mfem::Operator
   /// Named according to test variable.
   NamedFieldsMap<std::vector<std::shared_ptr<MFEMEssentialBC>>> _essential_bc_map;
 
+  // Operator handle for the jacobian matrix
   mutable mfem::OperatorHandle _jacobian;
-
   mfem::AssemblyLevel _assembly_level;
 
+  // Temporary vectors used for non-linear action
+  // assembly process
+  mutable mfem::BlockVector _trueBlockSol, _blockForces, _blockResidual;
+  Moose::MFEM::GridFunctions * _gfuncs;
+  mfem::Array<int> * _block_true_offsets = NULL;
+  mfem::Array<int> empty_tdof;
+  bool _non_linear = false;
+
 private:
   friend class EquationSystemProblemOperator;
   /// Disallowed inherited method
@@ -246,6 +269,29 @@ EquationSystem::ApplyDomainLFIntegrators(
   }
 }
 
+inline void
+EquationSystem::ApplyDomainNLAFIntegrators(
+    const std::string & test_var_name,
+    std::shared_ptr<mfem::ParLinearForm> form,
+    NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map)
+{
+  if (kernels_map.Has(test_var_name) && kernels_map.Get(test_var_name)->Has(test_var_name))
+  {
+    auto kernels = kernels_map.GetRef(test_var_name).GetRef(test_var_name);
+    for (auto & kernel : kernels)
+    {
+      mfem::LinearFormIntegrator * integ = kernel->createNLAIntegrator();
+      if (integ)
+      {
+        _non_linear = true;
+        kernel->isSubdomainRestricted()
+            ? form->AddDomainIntegrator(std::move(integ), kernel->getSubdomainMarkers())
+            : form->AddDomainIntegrator(std::move(integ));
+      }
+    }
+  }
+}
+
 template <class FormType>
 void
 EquationSystem::ApplyBoundaryBLFIntegrators(
@@ -301,6 +347,31 @@ EquationSystem::ApplyBoundaryLFIntegrators(
   }
 }
 
+inline void
+EquationSystem::ApplyBoundaryNLAFIntegrators(
+    const std::string & test_var_name,
+    std::shared_ptr<mfem::ParLinearForm> form,
+    NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
+        integrated_bc_map)
+{
+  if (integrated_bc_map.Has(test_var_name) &&
+      integrated_bc_map.Get(test_var_name)->Has(test_var_name))
+  {
+    auto bcs = integrated_bc_map.GetRef(test_var_name).GetRef(test_var_name);
+    for (auto & bc : bcs)
+    {
+      mfem::LinearFormIntegrator * integ = bc->createNLAIntegrator();
+      if (integ)
+      {
+        _non_linear = true;
+        bc->isBoundaryRestricted()
+            ? form->AddBoundaryIntegrator(std::move(integ), bc->getBoundaryMarkers())
+            : form->AddBoundaryIntegrator(std::move(integ));
+      }
+    }
+  }
+}
+
 } // namespace Moose::MFEM
 
 #endif
diff --git a/framework/include/mfem/executioners/MFEMProblemSolve.h b/framework/include/mfem/executioners/MFEMProblemSolve.h
index f41a692de044..606b35f57597 100644
--- a/framework/include/mfem/executioners/MFEMProblemSolve.h
+++ b/framework/include/mfem/executioners/MFEMProblemSolve.h
@@ -36,6 +36,10 @@ class MFEMProblemSolve : public SolveObject
 protected:
   MFEMProblem & _mfem_problem;
   std::vector<std::shared_ptr<Moose::MFEM::ProblemOperatorBase>> & _problem_operators;
+  unsigned int _nl_max_its;
+  mfem::real_t _nl_abs_tol;
+  mfem::real_t _nl_rel_tol;
+  unsigned int _print_level;
 };
 
 #endif
diff --git a/framework/include/mfem/kernels/MFEMKernel.h b/framework/include/mfem/kernels/MFEMKernel.h
index beeba72bf66a..a7db44b50a57 100644
--- a/framework/include/mfem/kernels/MFEMKernel.h
+++ b/framework/include/mfem/kernels/MFEMKernel.h
@@ -30,6 +30,7 @@ class MFEMKernel : public MFEMGeneralUserObject, public MFEMBlockRestrictable
   /// Create a new MFEM integrator to apply to the weak form. Ownership managed by the caller.
   virtual mfem::LinearFormIntegrator * createLFIntegrator() { return nullptr; }
   virtual mfem::BilinearFormIntegrator * createBFIntegrator() { return nullptr; }
+  virtual mfem::LinearFormIntegrator * createNLAIntegrator() { return nullptr; }
 
   /// Get name of the test variable labelling the weak form this kernel is added to
   const VariableName & getTestVariableName() const { return _test_var_name; }
diff --git a/framework/include/mfem/kernels/MFEMNLDiffusionKernel.h b/framework/include/mfem/kernels/MFEMNLDiffusionKernel.h
new file mode 100644
index 000000000000..4761367326d6
--- /dev/null
+++ b/framework/include/mfem/kernels/MFEMNLDiffusionKernel.h
@@ -0,0 +1,38 @@
+//* This file is part of the MOOSE framework
+//* https://mooseframework.inl.gov
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#ifdef MOOSE_MFEM_ENABLED
+
+#pragma once
+#include "MFEMKernel.h"
+
+/**
+ * \f[
+ * (k(u) \vec \nabla u, \vec \nabla v)
+ * \f]
+ */
+class MFEMNLDiffusionKernel : public MFEMKernel
+{
+public:
+  static InputParameters validParams();
+
+  MFEMNLDiffusionKernel(const InputParameters & parameters);
+
+  virtual mfem::BilinearFormIntegrator * createBFIntegrator() override;
+  virtual mfem::LinearFormIntegrator * createNLAIntegrator() override;
+
+protected:
+  mfem::Coefficient & _coef;
+  mfem::ScalarVectorProductCoefficient * _product_coef_res;
+  mfem::ScalarVectorProductCoefficient * _product_coef_jac;
+  mfem::SumIntegrator * _sum;
+  mfem::ConstantCoefficient * _one;
+};
+
+#endif
diff --git a/framework/include/mfem/problem/MFEMProblem.h b/framework/include/mfem/problem/MFEMProblem.h
index ce12509f2276..3b8346242a79 100644
--- a/framework/include/mfem/problem/MFEMProblem.h
+++ b/framework/include/mfem/problem/MFEMProblem.h
@@ -5,7 +5,7 @@
 //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
 //*
 //* Licensed under LGPL 2.1, please see LICENSE for details
-//* https://www.gnu.org/licenses/lgpl-2.1.html
+//* https://www.gnu.org/licenses/lgpl-2.1.html*/
 
 #ifdef MOOSE_MFEM_ENABLED
 
@@ -192,7 +192,10 @@ class MFEMProblem : public ExternalProblem
    * Add the nonlinear solver to the system. TODO: allow user to specify solver options,
    * similar to the linear solvers.
    */
-  void addMFEMNonlinearSolver();
+  void addMFEMNonlinearSolver(unsigned int nl_max_its,
+                              mfem::real_t nl_abs_tol,
+                              mfem::real_t nl_rel_tol,
+                              unsigned int print_level);
 
   /**
    * Method used to get an mfem FEC depending on the variable family specified in the input file.
diff --git a/framework/src/mfem/equation_systems/ComplexEquationSystem.C b/framework/src/mfem/equation_systems/ComplexEquationSystem.C
index c8f4d4e894cf..9ac7bc4ac235 100644
--- a/framework/src/mfem/equation_systems/ComplexEquationSystem.C
+++ b/framework/src/mfem/equation_systems/ComplexEquationSystem.C
@@ -269,6 +269,15 @@ ComplexEquationSystem::FormSystemMatrix(mfem::OperatorHandle & op,
   op.Reset(mfem::HypreParMatrixFromBlocks(_h_blocks));
 }
 
+// Equation system Mult
+void
+ComplexEquationSystem::Mult(const mfem::Vector & x, mfem::Vector & residual) const
+{
+  _jacobian->Mult(x, residual);
+  x.HostRead();
+  residual.HostRead();
+}
+
 void
 ComplexEquationSystem::RecoverComplexFEMSolution(
     mfem::BlockVector & trueX,
diff --git a/framework/src/mfem/equation_systems/EquationSystem.C b/framework/src/mfem/equation_systems/EquationSystem.C
index ec95320f8ed2..487d71df9405 100644
--- a/framework/src/mfem/equation_systems/EquationSystem.C
+++ b/framework/src/mfem/equation_systems/EquationSystem.C
@@ -5,7 +5,7 @@
 //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
 //*
 //* Licensed under LGPL 2.1, please see LICENSE for details
-//* https://www.gnu.org/licenses/lgpl-2.1.html
+//* https://www.gnu.org/licenses/lgpl-2.1.html */
 
 #ifdef MOOSE_MFEM_ENABLED
 
@@ -15,7 +15,14 @@
 namespace Moose::MFEM
 {
 
-EquationSystem::~EquationSystem() { DeleteAllBlocks(); }
+EquationSystem::~EquationSystem()
+{
+  if (_gfuncs != NULL)
+    delete _gfuncs;
+  if (_block_true_offsets != NULL)
+    delete _block_true_offsets;
+  DeleteAllBlocks();
+}
 
 void
 EquationSystem::DeleteAllBlocks()
@@ -159,6 +166,22 @@ EquationSystem::Init(Moose::MFEM::GridFunctions & gridfunctions,
   for (auto & eliminated_var_name : _eliminated_var_names)
     _eliminated_variables.Register(eliminated_var_name,
                                    gridfunctions.GetShared(eliminated_var_name));
+
+  // Get a reference to the GridFunctions
+  _gfuncs = new Moose::MFEM::GridFunctions(gridfunctions);
+
+  // Build the temporary BlockVectors
+  _block_true_offsets = new mfem::Array<int>(_trial_var_names.size() + 1);
+  (*_block_true_offsets)[0] = 0;
+
+  for (unsigned I = 0; I < _trial_var_names.size(); I++)
+  {
+    (*_block_true_offsets)[I + 1] = _gfuncs->Get(_trial_var_names.at(I))->ParFESpace()->TrueVSize();
+  }
+  _block_true_offsets->PartialSum();
+  _trueBlockSol.Update(*_block_true_offsets);
+  _blockForces.Update(*_block_true_offsets);
+  _blockResidual.Update(*_block_true_offsets);
 }
 
 void
@@ -322,6 +345,35 @@ EquationSystem::FormSystemMatrix(mfem::OperatorHandle & op,
   op.Reset(mfem::HypreParMatrixFromBlocks(_h_blocks));
 }
 
+void
+EquationSystem::ReassembleJacobian(mfem::BlockVector & x, mfem::BlockVector & rhs)
+{
+  // Reassemble all the Forms
+  for (const auto I : index_range(_test_var_names))
+  {
+    auto test_var_name = _test_var_names.at(I);
+    _blfs.GetShared(test_var_name)->Update();
+    _blfs.GetShared(test_var_name)->Assemble();
+    if (_mblfs.Has(test_var_name))
+    {
+      for (const auto J : index_range(_coupled_var_names))
+      {
+        auto coupled_var_name = _coupled_var_names.at(J);
+        if (_mblfs.Get(test_var_name)->Has(coupled_var_name))
+        {
+          _mblfs.GetShared(test_var_name)->GetShared(coupled_var_name)->Update();
+          _mblfs.GetShared(test_var_name)->GetShared(coupled_var_name)->Assemble();
+        }
+      }
+    }
+  }
+
+  // Form the system matrix
+  // This uses dummy arguments
+  // for the vectors
+  FormLinearSystem(_jacobian, x, rhs);
+}
+
 void
 EquationSystem::BuildJacobian(mfem::BlockVector & trueX, mfem::BlockVector & trueRHS)
 {
@@ -331,10 +383,89 @@ EquationSystem::BuildJacobian(mfem::BlockVector & trueX, mfem::BlockVector & tru
 }
 
 void
-EquationSystem::Mult(const mfem::Vector & x, mfem::Vector & residual) const
+EquationSystem::UpdateJacobian() const
 {
-  _jacobian->Mult(x, residual);
-  x.HostRead();
+
+  for (unsigned int i = 0; i < _test_var_names.size(); i++)
+  {
+    auto & test_var_name = _test_var_names.at(i);
+    auto blf = _blfs.Get(test_var_name);
+    blf->Update();
+    blf->Assemble();
+  }
+
+  // Form off-diagonal blocks
+  for (unsigned int i = 0; i < _test_var_names.size(); i++)
+  {
+    auto test_var_name = _test_var_names.at(i);
+    for (unsigned int j = 0; j < _test_var_names.size(); j++)
+    {
+      auto trial_var_name = _test_var_names.at(j);
+      if (_mblfs.Has(test_var_name) && _mblfs.Get(test_var_name)->Has(trial_var_name))
+      {
+        auto mblf = _mblfs.Get(test_var_name)->Get(trial_var_name);
+        mblf->Update();
+        mblf->Assemble();
+      }
+    }
+  }
+}
+
+void
+CopyVec(const mfem::Vector & x, mfem::Vector & y)
+{
+  y = x;
+}
+
+void
+EquationSystem::Mult(const mfem::Vector & sol, mfem::Vector & residual) const
+{
+  static_cast<mfem::Vector &>(_trueBlockSol) = sol;
+  for (unsigned int i = 0; i < _trial_var_names.size(); i++)
+  {
+    auto & trial_var_name = _trial_var_names.at(i);
+    _trueBlockSol.GetBlock(i).SyncAliasMemory(_trueBlockSol);
+    _gfuncs->Get(trial_var_name)->Distribute(&(_trueBlockSol.GetBlock(i)));
+  }
+
+  if (_non_linear)
+  {
+    _blockResidual = 0.0;
+    UpdateJacobian();
+
+    for (unsigned int i = 0; i < _test_var_names.size(); i++)
+    {
+      auto & test_var_name = _test_var_names.at(i);
+      int offset = _blockResidual.GetBlock(i).Size();
+      mfem::Vector b(offset);
+
+      auto lf = _lfs.GetShared(test_var_name);
+      lf->Assemble();
+      lf->ParallelAssemble(b);
+      b.SyncAliasMemory(b);
+
+      auto nlf = _nlAs.GetShared(test_var_name);
+      nlf->Assemble();
+      nlf->ParallelAssemble(_blockResidual.GetBlock(i));
+
+      _blockResidual.GetBlock(i) -= b;
+      _blockResidual.GetBlock(i) *= -1;
+
+      _blockResidual.GetBlock(i).SetSubVector(_ess_tdof_lists.at(i), 0.0);
+      _blockResidual.GetBlock(i).SyncAliasMemory(_blockResidual);
+    }
+
+    residual = static_cast<mfem::Vector &>(_blockResidual);
+    const_cast<EquationSystem *>(this)->FormLinearSystem(_jacobian, _trueBlockSol, _blockResidual);
+    residual *= -1.0;
+  }
+  else
+  {
+    residual = 0.0;
+    _jacobian->Mult(sol, residual);
+  }
+
+  sol.HostRead();
   residual.HostRead();
 }
 
@@ -384,6 +515,26 @@ EquationSystem::BuildLinearForms()
   EliminateCoupledVariables();
 }
 
+void
+EquationSystem::BuildNonLinearActionForms()
+{
+  // Register non-linear Action forms
+  for (const auto i : index_range(_test_var_names))
+  {
+    auto test_var_name = _test_var_names.at(i);
+    _nlAs.Register(test_var_name, std::make_shared<mfem::ParLinearForm>(_test_pfespaces.at(i)));
+    _nlAs.GetRef(test_var_name) = 0.0;
+  }
+
+  for (auto & test_var_name : _test_var_names)
+  {
+    // Apply kernels
+    auto nlA = _nlAs.GetShared(test_var_name);
+    ApplyDomainNLAFIntegrators(test_var_name, nlA, _kernels_map);
+    ApplyBoundaryNLAFIntegrators(test_var_name, nlA, _integrated_bc_map);
+  }
+}
+
 void
 EquationSystem::BuildBilinearForms()
 {
@@ -452,6 +603,7 @@ EquationSystem::BuildEquationSystem()
   BuildBilinearForms();
   BuildMixedBilinearForms();
   BuildLinearForms();
+  BuildNonLinearActionForms();
 }
 
 } // namespace Moose::MFEM
diff --git a/framework/src/mfem/executioners/MFEMProblemSolve.C b/framework/src/mfem/executioners/MFEMProblemSolve.C
index 7e54771df137..a183cafb8dd2 100644
--- a/framework/src/mfem/executioners/MFEMProblemSolve.C
+++ b/framework/src/mfem/executioners/MFEMProblemSolve.C
@@ -18,6 +18,10 @@ MFEMProblemSolve::validParams()
 {
   InputParameters params = emptyInputParameters();
   params.addClassDescription("Solve object for MFEM problems.");
+  params.addParam<unsigned int>("nl_max_its", 1, "Max Nonlinear Iterations");
+  params.addParam<Real>("nl_abs_tol", 1.0e-50, "Nonlinear Absolute Tolerance");
+  params.addParam<Real>("nl_rel_tol", 1.0e-8, "Nonlinear Relative Tolerance");
+  params.addParam<unsigned int>("print_level", 1, "Print level");
   params.addParam<std::string>("device", "Run app on the chosen device.");
   MooseEnum assembly_levels("legacy full element partial none", "legacy", true);
   params.addParam<MooseEnum>("assembly_level", assembly_levels, "Matrix assembly level.");
@@ -29,7 +33,11 @@ MFEMProblemSolve::MFEMProblemSolve(
     std::vector<std::shared_ptr<Moose::MFEM::ProblemOperatorBase>> & problem_operators)
   : SolveObject(ex),
     _mfem_problem(dynamic_cast<MFEMProblem &>(_problem)),
-    _problem_operators(problem_operators)
+    _problem_operators(problem_operators),
+    _nl_max_its(getParam<unsigned int>("nl_max_its")),
+    _nl_abs_tol(getParam<mfem::real_t>("nl_abs_tol")),
+    _nl_rel_tol(getParam<mfem::real_t>("nl_rel_tol")),
+    _print_level(getParam<unsigned int>("print_level"))
 {
   if (const auto compute_device = _app.getComputeDevice())
     _app.setMFEMDevice(*compute_device, Moose::PassKey<MFEMProblemSolve>());
@@ -38,6 +46,7 @@ MFEMProblemSolve::MFEMProblemSolve(
                        : _app.isUltimateMaster() ? "cpu"
                                                  : "",
                        Moose::PassKey<MFEMProblemSolve>());
+  _mfem_problem.addMFEMNonlinearSolver(_nl_max_its, _nl_abs_tol, _nl_rel_tol, _print_level);
 }
 
 bool
diff --git a/framework/src/mfem/kernels/MFEMNLDiffusionKernel.C b/framework/src/mfem/kernels/MFEMNLDiffusionKernel.C
new file mode 100644
index 000000000000..e5417ac66900
--- /dev/null
+++ b/framework/src/mfem/kernels/MFEMNLDiffusionKernel.C
@@ -0,0 +1,65 @@
+//* This file is part of the MOOSE framework
+//* https://mooseframework.inl.gov
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#ifdef MOOSE_MFEM_ENABLED
+
+#include "MFEMNLDiffusionKernel.h"
+#include "MFEMProblem.h"
+
+registerMooseObject("MooseApp", MFEMNLDiffusionKernel);
+
+InputParameters
+MFEMNLDiffusionKernel::validParams()
+{
+  InputParameters params = MFEMKernel::validParams();
+  params.addClassDescription("Adds the domain integrator for integrating the non-linear action"
+                             "$(k(u)\\vec\\nabla v, \\vec\\nabla v)_\\Omega$"
+                             "Adds the domain integrator to an MFEM problem for the bilinear form "
+                             "$(k(u)\\vec\\nable v, \\vec\\nabla v)_\\Omega + (k'(u) v, "
+                             "\\vec\\nabla u \\vec\\nabla v)_\\Omega$ "
+                             "The above terms arises from the weak form of the non-linear operator "
+                             "$- \\vec\\nabla \\cdot ( k(u) \\vec\\nabla u)$.");
+  params.addParam<MFEMScalarCoefficientName>(
+      "coefficient", "1.", "Name of property for MixedScalarWeakDivergence coefficient k.");
+  return params;
+}
+
+MFEMNLDiffusionKernel::MFEMNLDiffusionKernel(const InputParameters & parameters)
+  : MFEMKernel(parameters), _coef(getScalarCoefficient("coefficient"))
+// FIXME: The MFEM bilinear form can also handle vector and matrix
+// coefficients, so ideally we'd handle all three too.
+{
+  // declares GradientGridFunctionCoefficient
+  getMFEMProblem().getCoefficients().declareVector<mfem::GradientGridFunctionCoefficient>(
+      name(), getMFEMProblem().getProblemData().gridfunctions.Get(_test_var_name));
+}
+
+mfem::BilinearFormIntegrator *
+MFEMNLDiffusionKernel::createBFIntegrator()
+{
+  _sum = new mfem::SumIntegrator;
+  _sum->AddIntegrator(new mfem::DiffusionIntegrator(_coef));
+  mfem::VectorCoefficient & vec_coef =
+      getMFEMProblem().getCoefficients().getVectorCoefficient(name());
+  _one = new mfem::ConstantCoefficient(1.0);
+  _product_coef_jac = new mfem::ScalarVectorProductCoefficient(*_one, vec_coef);
+  _sum->AddIntegrator(new mfem::MixedScalarWeakDivergenceIntegrator(*_product_coef_jac));
+  return _sum;
+}
+
+mfem::LinearFormIntegrator *
+MFEMNLDiffusionKernel::createNLAIntegrator()
+{
+  mfem::VectorCoefficient & vec_coef =
+      getMFEMProblem().getCoefficients().getVectorCoefficient(name());
+  _product_coef_res = new mfem::ScalarVectorProductCoefficient(_coef, vec_coef);
+  return new mfem::DomainLFGradIntegrator(*_product_coef_res);
+}
+
+#endif
diff --git a/framework/src/mfem/problem/MFEMProblem.C b/framework/src/mfem/problem/MFEMProblem.C
index e45150feffe0..f059953091ad 100644
--- a/framework/src/mfem/problem/MFEMProblem.C
+++ b/framework/src/mfem/problem/MFEMProblem.C
@@ -48,7 +48,6 @@ void
 MFEMProblem::initialSetup()
 {
   FEProblemBase::initialSetup();
-  addMFEMNonlinearSolver();
 }
 
 void
@@ -81,15 +80,18 @@ MFEMProblem::addMFEMSolver(const std::string & user_object_name,
 }
 
 void
-MFEMProblem::addMFEMNonlinearSolver()
+MFEMProblem::addMFEMNonlinearSolver(unsigned int nl_max_its,
+                                    mfem::real_t nl_abs_tol,
+                                    mfem::real_t nl_rel_tol,
+                                    unsigned int print_level)
 {
   auto nl_solver = std::make_shared<mfem::NewtonSolver>(getComm());
 
   // Defaults to one iteration, without further nonlinear iterations
-  nl_solver->SetRelTol(0.0);
-  nl_solver->SetAbsTol(0.0);
-  nl_solver->SetMaxIter(1);
-
+  nl_solver->SetRelTol(nl_rel_tol);
+  nl_solver->SetAbsTol(nl_abs_tol);
+  nl_solver->SetMaxIter(nl_max_its);
+  nl_solver->SetPrintLevel(print_level);
   getProblemData().nonlinear_solver = nl_solver;
 }
 
diff --git a/framework/src/mfem/problem_operators/EquationSystemProblemOperator.C b/framework/src/mfem/problem_operators/EquationSystemProblemOperator.C
index fd7358fe3ec9..1254b72c829f 100644
--- a/framework/src/mfem/problem_operators/EquationSystemProblemOperator.C
+++ b/framework/src/mfem/problem_operators/EquationSystemProblemOperator.C
@@ -30,17 +30,27 @@ EquationSystemProblemOperator::Init(mfem::BlockVector & X)
 void
 EquationSystemProblemOperator::Solve()
 {
-  GetEquationSystem()->BuildJacobian(_true_x, _true_rhs);
 
+  GetEquationSystem()->BuildJacobian(_true_x, _true_rhs);
   if (_problem_data.jacobian_solver->isLOR() && GetEquationSystem()->GetTestVarNames().size() > 1)
     mooseError("LOR solve is only supported for single-variable systems");
   _problem_data.jacobian_solver->updateSolver(
       *GetEquationSystem()->_blfs.Get(GetEquationSystem()->GetTestVarNames().at(0)),
       GetEquationSystem()->_ess_tdof_lists.at(0));
-
+  /*
+    _problem_data.nonlinear_solver->SetSolver(_problem_data.jacobian_solver->getSolver());
+    _problem_data.nonlinear_solver->SetOperator(*GetEquationSystem());
+    _problem_data.nonlinear_solver->Mult(_true_rhs, _true_x);
+  */
+  mfem::Vector zero_vec(_true_rhs.Size());
+  zero_vec = 0.0;
   _problem_data.nonlinear_solver->SetSolver(_problem_data.jacobian_solver->getSolver());
   _problem_data.nonlinear_solver->SetOperator(*GetEquationSystem());
-  _problem_data.nonlinear_solver->Mult(_true_rhs, _true_x);
+
+  if (!(GetEquationSystem()->_non_linear))
+    _problem_data.nonlinear_solver->Mult(_true_rhs, _true_x);
+  else
+    _problem_data.nonlinear_solver->Mult(zero_vec, _true_x);
 
   GetEquationSystem()->RecoverFEMSolution(
       _true_x, _problem_data.gridfunctions, _problem_data.cmplx_gridfunctions);
diff --git a/test/tests/mfem/kernels/gold/OutputData/NLDiffusion/Run0/Cycle000000/data.pvtu b/test/tests/mfem/kernels/gold/OutputData/NLDiffusion/Run0/Cycle000000/data.pvtu
new file mode 100644
index 000000000000..8efaaa4e6a2e
--- /dev/null
+++ b/test/tests/mfem/kernels/gold/OutputData/NLDiffusion/Run0/Cycle000000/data.pvtu
@@ -0,0 +1,20 @@
+<?xml version="1.0"?>
+<VTKFile type="PUnstructuredGrid" version ="2.2" byte_order="LittleEndian">
+<PUnstructuredGrid GhostLevel="0">
+<PPoints>
+	<PDataArray type="Float64"  Name="Points" NumberOfComponents="3" format="ascii"/>
+</PPoints>
+<PCells>
+	<PDataArray type="Int32"  Name="connectivity" NumberOfComponents="1" format="ascii"/>
+	<PDataArray type="Int32"  Name="offsets"      NumberOfComponents="1" format="ascii"/>
+	<PDataArray type="UInt8"  Name="types"        NumberOfComponents="1" format="ascii"/>
+</PCells>
+<PPointData>
+<PDataArray type="Float64" Name="concentration" NumberOfComponents="1"  format="ascii" />
+</PPointData>
+<PCellData>
+	<PDataArray type="Int32" Name="attribute" NumberOfComponents="1" format="ascii"/>
+</PCellData>
+<Piece Source="proc000000.vtu"/>
+</PUnstructuredGrid>
+</VTKFile>
diff --git a/test/tests/mfem/kernels/gold/OutputData/NLDiffusion/Run0/Cycle000000/proc000000.vtu b/test/tests/mfem/kernels/gold/OutputData/NLDiffusion/Run0/Cycle000000/proc000000.vtu
new file mode 100644
index 000000000000..449315f29594
--- /dev/null
+++ b/test/tests/mfem/kernels/gold/OutputData/NLDiffusion/Run0/Cycle000000/proc000000.vtu
@@ -0,0 +1,1226 @@
+<VTKFile type="UnstructuredGrid" version="2.2" byte_order="LittleEndian">
+<UnstructuredGrid>
+<Piece NumberOfPoints="400" NumberOfCells="100">
+<Points>
+<DataArray type="Float64" NumberOfComponents="3" format="ascii">
+0 0 0
+0.1 0 0
+0 0.1 0
+0.1 0.1 0
+0.1 0 0
+0.2 0 0
+0.1 0.1 0
+0.2 0.1 0
+0.2 0 0
+0.3 0 0
+0.2 0.1 0
+0.3 0.1 0
+0.3 0 0
+0.4 0 0
+0.3 0.1 0
+0.4 0.1 0
+0.4 0 0
+0.5 0 0
+0.4 0.1 0
+0.5 0.1 0
+0.5 0 0
+0.6 0 0
+0.5 0.1 0
+0.6 0.1 0
+0.6 0 0
+0.7 0 0
+0.6 0.1 0
+0.7 0.1 0
+0.7 0 0
+0.8 0 0
+0.7 0.1 0
+0.8 0.1 0
+0.8 0 0
+0.9 0 0
+0.8 0.1 0
+0.9 0.1 0
+0.9 0 0
+1 0 0
+0.9 0.1 0
+1 0.1 0
+0 0.1 0
+0.1 0.1 0
+0 0.2 0
+0.1 0.2 0
+0.1 0.1 0
+0.2 0.1 0
+0.1 0.2 0
+0.2 0.2 0
+0.2 0.1 0
+0.3 0.1 0
+0.2 0.2 0
+0.3 0.2 0
+0.3 0.1 0
+0.4 0.1 0
+0.3 0.2 0
+0.4 0.2 0
+0.4 0.1 0
+0.5 0.1 0
+0.4 0.2 0
+0.5 0.2 0
+0.5 0.1 0
+0.6 0.1 0
+0.5 0.2 0
+0.6 0.2 0
+0.6 0.1 0
+0.7 0.1 0
+0.6 0.2 0
+0.7 0.2 0
+0.7 0.1 0
+0.8 0.1 0
+0.7 0.2 0
+0.8 0.2 0
+0.8 0.1 0
+0.9 0.1 0
+0.8 0.2 0
+0.9 0.2 0
+0.9 0.1 0
+1 0.1 0
+0.9 0.2 0
+1 0.2 0
+0 0.2 0
+0.1 0.2 0
+0 0.3 0
+0.1 0.3 0
+0.1 0.2 0
+0.2 0.2 0
+0.1 0.3 0
+0.2 0.3 0
+0.2 0.2 0
+0.3 0.2 0
+0.2 0.3 0
+0.3 0.3 0
+0.3 0.2 0
+0.4 0.2 0
+0.3 0.3 0
+0.4 0.3 0
+0.4 0.2 0
+0.5 0.2 0
+0.4 0.3 0
+0.5 0.3 0
+0.5 0.2 0
+0.6 0.2 0
+0.5 0.3 0
+0.6 0.3 0
+0.6 0.2 0
+0.7 0.2 0
+0.6 0.3 0
+0.7 0.3 0
+0.7 0.2 0
+0.8 0.2 0
+0.7 0.3 0
+0.8 0.3 0
+0.8 0.2 0
+0.9 0.2 0
+0.8 0.3 0
+0.9 0.3 0
+0.9 0.2 0
+1 0.2 0
+0.9 0.3 0
+1 0.3 0
+0 0.3 0
+0.1 0.3 0
+0 0.4 0
+0.1 0.4 0
+0.1 0.3 0
+0.2 0.3 0
+0.1 0.4 0
+0.2 0.4 0
+0.2 0.3 0
+0.3 0.3 0
+0.2 0.4 0
+0.3 0.4 0
+0.3 0.3 0
+0.4 0.3 0
+0.3 0.4 0
+0.4 0.4 0
+0.4 0.3 0
+0.5 0.3 0
+0.4 0.4 0
+0.5 0.4 0
+0.5 0.3 0
+0.6 0.3 0
+0.5 0.4 0
+0.6 0.4 0
+0.6 0.3 0
+0.7 0.3 0
+0.6 0.4 0
+0.7 0.4 0
+0.7 0.3 0
+0.8 0.3 0
+0.7 0.4 0
+0.8 0.4 0
+0.8 0.3 0
+0.9 0.3 0
+0.8 0.4 0
+0.9 0.4 0
+0.9 0.3 0
+1 0.3 0
+0.9 0.4 0
+1 0.4 0
+0 0.4 0
+0.1 0.4 0
+0 0.5 0
+0.1 0.5 0
+0.1 0.4 0
+0.2 0.4 0
+0.1 0.5 0
+0.2 0.5 0
+0.2 0.4 0
+0.3 0.4 0
+0.2 0.5 0
+0.3 0.5 0
+0.3 0.4 0
+0.4 0.4 0
+0.3 0.5 0
+0.4 0.5 0
+0.4 0.4 0
+0.5 0.4 0
+0.4 0.5 0
+0.5 0.5 0
+0.5 0.4 0
+0.6 0.4 0
+0.5 0.5 0
+0.6 0.5 0
+0.6 0.4 0
+0.7 0.4 0
+0.6 0.5 0
+0.7 0.5 0
+0.7 0.4 0
+0.8 0.4 0
+0.7 0.5 0
+0.8 0.5 0
+0.8 0.4 0
+0.9 0.4 0
+0.8 0.5 0
+0.9 0.5 0
+0.9 0.4 0
+1 0.4 0
+0.9 0.5 0
+1 0.5 0
+0 0.5 0
+0.1 0.5 0
+0 0.6 0
+0.1 0.6 0
+0.1 0.5 0
+0.2 0.5 0
+0.1 0.6 0
+0.2 0.6 0
+0.2 0.5 0
+0.3 0.5 0
+0.2 0.6 0
+0.3 0.6 0
+0.3 0.5 0
+0.4 0.5 0
+0.3 0.6 0
+0.4 0.6 0
+0.4 0.5 0
+0.5 0.5 0
+0.4 0.6 0
+0.5 0.6 0
+0.5 0.5 0
+0.6 0.5 0
+0.5 0.6 0
+0.6 0.6 0
+0.6 0.5 0
+0.7 0.5 0
+0.6 0.6 0
+0.7 0.6 0
+0.7 0.5 0
+0.8 0.5 0
+0.7 0.6 0
+0.8 0.6 0
+0.8 0.5 0
+0.9 0.5 0
+0.8 0.6 0
+0.9 0.6 0
+0.9 0.5 0
+1 0.5 0
+0.9 0.6 0
+1 0.6 0
+0 0.6 0
+0.1 0.6 0
+0 0.7 0
+0.1 0.7 0
+0.1 0.6 0
+0.2 0.6 0
+0.1 0.7 0
+0.2 0.7 0
+0.2 0.6 0
+0.3 0.6 0
+0.2 0.7 0
+0.3 0.7 0
+0.3 0.6 0
+0.4 0.6 0
+0.3 0.7 0
+0.4 0.7 0
+0.4 0.6 0
+0.5 0.6 0
+0.4 0.7 0
+0.5 0.7 0
+0.5 0.6 0
+0.6 0.6 0
+0.5 0.7 0
+0.6 0.7 0
+0.6 0.6 0
+0.7 0.6 0
+0.6 0.7 0
+0.7 0.7 0
+0.7 0.6 0
+0.8 0.6 0
+0.7 0.7 0
+0.8 0.7 0
+0.8 0.6 0
+0.9 0.6 0
+0.8 0.7 0
+0.9 0.7 0
+0.9 0.6 0
+1 0.6 0
+0.9 0.7 0
+1 0.7 0
+0 0.7 0
+0.1 0.7 0
+0 0.8 0
+0.1 0.8 0
+0.1 0.7 0
+0.2 0.7 0
+0.1 0.8 0
+0.2 0.8 0
+0.2 0.7 0
+0.3 0.7 0
+0.2 0.8 0
+0.3 0.8 0
+0.3 0.7 0
+0.4 0.7 0
+0.3 0.8 0
+0.4 0.8 0
+0.4 0.7 0
+0.5 0.7 0
+0.4 0.8 0
+0.5 0.8 0
+0.5 0.7 0
+0.6 0.7 0
+0.5 0.8 0
+0.6 0.8 0
+0.6 0.7 0
+0.7 0.7 0
+0.6 0.8 0
+0.7 0.8 0
+0.7 0.7 0
+0.8 0.7 0
+0.7 0.8 0
+0.8 0.8 0
+0.8 0.7 0
+0.9 0.7 0
+0.8 0.8 0
+0.9 0.8 0
+0.9 0.7 0
+1 0.7 0
+0.9 0.8 0
+1 0.8 0
+0 0.8 0
+0.1 0.8 0
+0 0.9 0
+0.1 0.9 0
+0.1 0.8 0
+0.2 0.8 0
+0.1 0.9 0
+0.2 0.9 0
+0.2 0.8 0
+0.3 0.8 0
+0.2 0.9 0
+0.3 0.9 0
+0.3 0.8 0
+0.4 0.8 0
+0.3 0.9 0
+0.4 0.9 0
+0.4 0.8 0
+0.5 0.8 0
+0.4 0.9 0
+0.5 0.9 0
+0.5 0.8 0
+0.6 0.8 0
+0.5 0.9 0
+0.6 0.9 0
+0.6 0.8 0
+0.7 0.8 0
+0.6 0.9 0
+0.7 0.9 0
+0.7 0.8 0
+0.8 0.8 0
+0.7 0.9 0
+0.8 0.9 0
+0.8 0.8 0
+0.9 0.8 0
+0.8 0.9 0
+0.9 0.9 0
+0.9 0.8 0
+1 0.8 0
+0.9 0.9 0
+1 0.9 0
+0 0.9 0
+0.1 0.9 0
+0 1 0
+0.1 1 0
+0.1 0.9 0
+0.2 0.9 0
+0.1 1 0
+0.2 1 0
+0.2 0.9 0
+0.3 0.9 0
+0.2 1 0
+0.3 1 0
+0.3 0.9 0
+0.4 0.9 0
+0.3 1 0
+0.4 1 0
+0.4 0.9 0
+0.5 0.9 0
+0.4 1 0
+0.5 1 0
+0.5 0.9 0
+0.6 0.9 0
+0.5 1 0
+0.6 1 0
+0.6 0.9 0
+0.7 0.9 0
+0.6 1 0
+0.7 1 0
+0.7 0.9 0
+0.8 0.9 0
+0.7 1 0
+0.8 1 0
+0.8 0.9 0
+0.9 0.9 0
+0.8 1 0
+0.9 1 0
+0.9 0.9 0
+1 0.9 0
+0.9 1 0
+1 1 0
+</DataArray>
+</Points>
+<Cells>
+<DataArray type="Int32" Name="connectivity" format="ascii">
+0 1 3 2 
+4 5 7 6 
+8 9 11 10 
+12 13 15 14 
+16 17 19 18 
+20 21 23 22 
+24 25 27 26 
+28 29 31 30 
+32 33 35 34 
+36 37 39 38 
+40 41 43 42 
+44 45 47 46 
+48 49 51 50 
+52 53 55 54 
+56 57 59 58 
+60 61 63 62 
+64 65 67 66 
+68 69 71 70 
+72 73 75 74 
+76 77 79 78 
+80 81 83 82 
+84 85 87 86 
+88 89 91 90 
+92 93 95 94 
+96 97 99 98 
+100 101 103 102 
+104 105 107 106 
+108 109 111 110 
+112 113 115 114 
+116 117 119 118 
+120 121 123 122 
+124 125 127 126 
+128 129 131 130 
+132 133 135 134 
+136 137 139 138 
+140 141 143 142 
+144 145 147 146 
+148 149 151 150 
+152 153 155 154 
+156 157 159 158 
+160 161 163 162 
+164 165 167 166 
+168 169 171 170 
+172 173 175 174 
+176 177 179 178 
+180 181 183 182 
+184 185 187 186 
+188 189 191 190 
+192 193 195 194 
+196 197 199 198 
+200 201 203 202 
+204 205 207 206 
+208 209 211 210 
+212 213 215 214 
+216 217 219 218 
+220 221 223 222 
+224 225 227 226 
+228 229 231 230 
+232 233 235 234 
+236 237 239 238 
+240 241 243 242 
+244 245 247 246 
+248 249 251 250 
+252 253 255 254 
+256 257 259 258 
+260 261 263 262 
+264 265 267 266 
+268 269 271 270 
+272 273 275 274 
+276 277 279 278 
+280 281 283 282 
+284 285 287 286 
+288 289 291 290 
+292 293 295 294 
+296 297 299 298 
+300 301 303 302 
+304 305 307 306 
+308 309 311 310 
+312 313 315 314 
+316 317 319 318 
+320 321 323 322 
+324 325 327 326 
+328 329 331 330 
+332 333 335 334 
+336 337 339 338 
+340 341 343 342 
+344 345 347 346 
+348 349 351 350 
+352 353 355 354 
+356 357 359 358 
+360 361 363 362 
+364 365 367 366 
+368 369 371 370 
+372 373 375 374 
+376 377 379 378 
+380 381 383 382 
+384 385 387 386 
+388 389 391 390 
+392 393 395 394 
+396 397 399 398 
+</DataArray>
+<DataArray type="Int32" Name="offsets" format="ascii">
+4
+8
+12
+16
+20
+24
+28
+32
+36
+40
+44
+48
+52
+56
+60
+64
+68
+72
+76
+80
+84
+88
+92
+96
+100
+104
+108
+112
+116
+120
+124
+128
+132
+136
+140
+144
+148
+152
+156
+160
+164
+168
+172
+176
+180
+184
+188
+192
+196
+200
+204
+208
+212
+216
+220
+224
+228
+232
+236
+240
+244
+248
+252
+256
+260
+264
+268
+272
+276
+280
+284
+288
+292
+296
+300
+304
+308
+312
+316
+320
+324
+328
+332
+336
+340
+344
+348
+352
+356
+360
+364
+368
+372
+376
+380
+384
+388
+392
+396
+400
+</DataArray>
+<DataArray type="UInt8" Name="types" format="ascii">
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+</DataArray>
+</Cells>
+<CellData Scalars="attribute">
+<DataArray type="Int32" Name="attribute" format="ascii">
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+</DataArray>
+</CellData>
+<PointData >
+<DataArray type="Float64" Name="concentration" NumberOfComponents="1"  format="ascii" >
+1
+1
+1.2
+1.2
+1
+1
+1.2
+1.2
+1
+1
+1.2
+1.2
+1
+1
+1.2
+1.2
+1
+1
+1.2
+1.2
+1
+1
+1.2
+1.2
+1
+1
+1.2
+1.2
+1
+1
+1.2
+1.2
+1
+1
+1.2
+1.2
+1
+1
+1.2
+1.2
+1.2
+1.2
+1.4
+1.4
+1.2
+1.2
+1.4
+1.4
+1.2
+1.2
+1.4
+1.4
+1.2
+1.2
+1.4
+1.4
+1.2
+1.2
+1.4
+1.4
+1.2
+1.2
+1.4
+1.4
+1.2
+1.2
+1.4
+1.4
+1.2
+1.2
+1.4
+1.4
+1.2
+1.2
+1.4
+1.4
+1.2
+1.2
+1.4
+1.4
+1.4
+1.4
+1.6
+1.6
+1.4
+1.4
+1.6
+1.6
+1.4
+1.4
+1.6
+1.6
+1.4
+1.4
+1.6
+1.6
+1.4
+1.4
+1.6
+1.6
+1.4
+1.4
+1.6
+1.6
+1.4
+1.4
+1.6
+1.6
+1.4
+1.4
+1.6
+1.6
+1.4
+1.4
+1.6
+1.6
+1.4
+1.4
+1.6
+1.6
+1.6
+1.6
+1.8
+1.8
+1.6
+1.6
+1.8
+1.8
+1.6
+1.6
+1.8
+1.8
+1.6
+1.6
+1.8
+1.8
+1.6
+1.6
+1.8
+1.8
+1.6
+1.6
+1.8
+1.8
+1.6
+1.6
+1.8
+1.8
+1.6
+1.6
+1.8
+1.8
+1.6
+1.6
+1.8
+1.8
+1.6
+1.6
+1.8
+1.8
+1.8
+1.8
+2
+2
+1.8
+1.8
+2
+2
+1.8
+1.8
+2
+2
+1.8
+1.8
+2
+2
+1.8
+1.8
+2
+2
+1.8
+1.8
+2
+2
+1.8
+1.8
+2
+2
+1.8
+1.8
+2
+2
+1.8
+1.8
+2
+2
+1.8
+1.8
+2
+2
+2
+2
+2.2
+2.2
+2
+2
+2.2
+2.2
+2
+2
+2.2
+2.2
+2
+2
+2.2
+2.2
+2
+2
+2.2
+2.2
+2
+2
+2.2
+2.2
+2
+2
+2.2
+2.2
+2
+2
+2.2
+2.2
+2
+2
+2.2
+2.2
+2
+2
+2.2
+2.2
+2.2
+2.2
+2.4
+2.4
+2.2
+2.2
+2.4
+2.4
+2.2
+2.2
+2.4
+2.4
+2.2
+2.2
+2.4
+2.4
+2.2
+2.2
+2.4
+2.4
+2.2
+2.2
+2.4
+2.4
+2.2
+2.2
+2.4
+2.4
+2.2
+2.2
+2.4
+2.4
+2.2
+2.2
+2.4
+2.4
+2.2
+2.2
+2.4
+2.4
+2.4
+2.4
+2.6
+2.6
+2.4
+2.4
+2.6
+2.6
+2.4
+2.4
+2.6
+2.6
+2.4
+2.4
+2.6
+2.6
+2.4
+2.4
+2.6
+2.6
+2.4
+2.4
+2.6
+2.6
+2.4
+2.4
+2.6
+2.6
+2.4
+2.4
+2.6
+2.6
+2.4
+2.4
+2.6
+2.6
+2.4
+2.4
+2.6
+2.6
+2.6
+2.6
+2.8
+2.8
+2.6
+2.6
+2.8
+2.8
+2.6
+2.6
+2.8
+2.8
+2.6
+2.6
+2.8
+2.8
+2.6
+2.6
+2.8
+2.8
+2.6
+2.6
+2.8
+2.8
+2.6
+2.6
+2.8
+2.8
+2.6
+2.6
+2.8
+2.8
+2.6
+2.6
+2.8
+2.8
+2.6
+2.6
+2.8
+2.8
+2.8
+2.8
+3
+3
+2.8
+2.8
+3
+3
+2.8
+2.8
+3
+3
+2.8
+2.8
+3
+3
+2.8
+2.8
+3
+3
+2.8
+2.8
+3
+3
+2.8
+2.8
+3
+3
+2.8
+2.8
+3
+3
+2.8
+2.8
+3
+3
+2.8
+2.8
+3
+3
+</DataArray>
+</PointData>
+</Piece>
+</UnstructuredGrid>
+</VTKFile>
diff --git a/test/tests/mfem/kernels/gold/OutputData/NLDiffusion/Run0/Cycle000001/data.pvtu b/test/tests/mfem/kernels/gold/OutputData/NLDiffusion/Run0/Cycle000001/data.pvtu
new file mode 100644
index 000000000000..8efaaa4e6a2e
--- /dev/null
+++ b/test/tests/mfem/kernels/gold/OutputData/NLDiffusion/Run0/Cycle000001/data.pvtu
@@ -0,0 +1,20 @@
+<?xml version="1.0"?>
+<VTKFile type="PUnstructuredGrid" version ="2.2" byte_order="LittleEndian">
+<PUnstructuredGrid GhostLevel="0">
+<PPoints>
+	<PDataArray type="Float64"  Name="Points" NumberOfComponents="3" format="ascii"/>
+</PPoints>
+<PCells>
+	<PDataArray type="Int32"  Name="connectivity" NumberOfComponents="1" format="ascii"/>
+	<PDataArray type="Int32"  Name="offsets"      NumberOfComponents="1" format="ascii"/>
+	<PDataArray type="UInt8"  Name="types"        NumberOfComponents="1" format="ascii"/>
+</PCells>
+<PPointData>
+<PDataArray type="Float64" Name="concentration" NumberOfComponents="1"  format="ascii" />
+</PPointData>
+<PCellData>
+	<PDataArray type="Int32" Name="attribute" NumberOfComponents="1" format="ascii"/>
+</PCellData>
+<Piece Source="proc000000.vtu"/>
+</PUnstructuredGrid>
+</VTKFile>
diff --git a/test/tests/mfem/kernels/gold/OutputData/NLDiffusion/Run0/Cycle000001/proc000000.vtu b/test/tests/mfem/kernels/gold/OutputData/NLDiffusion/Run0/Cycle000001/proc000000.vtu
new file mode 100644
index 000000000000..f0b61ea4ec3e
--- /dev/null
+++ b/test/tests/mfem/kernels/gold/OutputData/NLDiffusion/Run0/Cycle000001/proc000000.vtu
@@ -0,0 +1,1226 @@
+<VTKFile type="UnstructuredGrid" version="2.2" byte_order="LittleEndian">
+<UnstructuredGrid>
+<Piece NumberOfPoints="400" NumberOfCells="100">
+<Points>
+<DataArray type="Float64" NumberOfComponents="3" format="ascii">
+0 0 0
+0.1 0 0
+0 0.1 0
+0.1 0.1 0
+0.1 0 0
+0.2 0 0
+0.1 0.1 0
+0.2 0.1 0
+0.2 0 0
+0.3 0 0
+0.2 0.1 0
+0.3 0.1 0
+0.3 0 0
+0.4 0 0
+0.3 0.1 0
+0.4 0.1 0
+0.4 0 0
+0.5 0 0
+0.4 0.1 0
+0.5 0.1 0
+0.5 0 0
+0.6 0 0
+0.5 0.1 0
+0.6 0.1 0
+0.6 0 0
+0.7 0 0
+0.6 0.1 0
+0.7 0.1 0
+0.7 0 0
+0.8 0 0
+0.7 0.1 0
+0.8 0.1 0
+0.8 0 0
+0.9 0 0
+0.8 0.1 0
+0.9 0.1 0
+0.9 0 0
+1 0 0
+0.9 0.1 0
+1 0.1 0
+0 0.1 0
+0.1 0.1 0
+0 0.2 0
+0.1 0.2 0
+0.1 0.1 0
+0.2 0.1 0
+0.1 0.2 0
+0.2 0.2 0
+0.2 0.1 0
+0.3 0.1 0
+0.2 0.2 0
+0.3 0.2 0
+0.3 0.1 0
+0.4 0.1 0
+0.3 0.2 0
+0.4 0.2 0
+0.4 0.1 0
+0.5 0.1 0
+0.4 0.2 0
+0.5 0.2 0
+0.5 0.1 0
+0.6 0.1 0
+0.5 0.2 0
+0.6 0.2 0
+0.6 0.1 0
+0.7 0.1 0
+0.6 0.2 0
+0.7 0.2 0
+0.7 0.1 0
+0.8 0.1 0
+0.7 0.2 0
+0.8 0.2 0
+0.8 0.1 0
+0.9 0.1 0
+0.8 0.2 0
+0.9 0.2 0
+0.9 0.1 0
+1 0.1 0
+0.9 0.2 0
+1 0.2 0
+0 0.2 0
+0.1 0.2 0
+0 0.3 0
+0.1 0.3 0
+0.1 0.2 0
+0.2 0.2 0
+0.1 0.3 0
+0.2 0.3 0
+0.2 0.2 0
+0.3 0.2 0
+0.2 0.3 0
+0.3 0.3 0
+0.3 0.2 0
+0.4 0.2 0
+0.3 0.3 0
+0.4 0.3 0
+0.4 0.2 0
+0.5 0.2 0
+0.4 0.3 0
+0.5 0.3 0
+0.5 0.2 0
+0.6 0.2 0
+0.5 0.3 0
+0.6 0.3 0
+0.6 0.2 0
+0.7 0.2 0
+0.6 0.3 0
+0.7 0.3 0
+0.7 0.2 0
+0.8 0.2 0
+0.7 0.3 0
+0.8 0.3 0
+0.8 0.2 0
+0.9 0.2 0
+0.8 0.3 0
+0.9 0.3 0
+0.9 0.2 0
+1 0.2 0
+0.9 0.3 0
+1 0.3 0
+0 0.3 0
+0.1 0.3 0
+0 0.4 0
+0.1 0.4 0
+0.1 0.3 0
+0.2 0.3 0
+0.1 0.4 0
+0.2 0.4 0
+0.2 0.3 0
+0.3 0.3 0
+0.2 0.4 0
+0.3 0.4 0
+0.3 0.3 0
+0.4 0.3 0
+0.3 0.4 0
+0.4 0.4 0
+0.4 0.3 0
+0.5 0.3 0
+0.4 0.4 0
+0.5 0.4 0
+0.5 0.3 0
+0.6 0.3 0
+0.5 0.4 0
+0.6 0.4 0
+0.6 0.3 0
+0.7 0.3 0
+0.6 0.4 0
+0.7 0.4 0
+0.7 0.3 0
+0.8 0.3 0
+0.7 0.4 0
+0.8 0.4 0
+0.8 0.3 0
+0.9 0.3 0
+0.8 0.4 0
+0.9 0.4 0
+0.9 0.3 0
+1 0.3 0
+0.9 0.4 0
+1 0.4 0
+0 0.4 0
+0.1 0.4 0
+0 0.5 0
+0.1 0.5 0
+0.1 0.4 0
+0.2 0.4 0
+0.1 0.5 0
+0.2 0.5 0
+0.2 0.4 0
+0.3 0.4 0
+0.2 0.5 0
+0.3 0.5 0
+0.3 0.4 0
+0.4 0.4 0
+0.3 0.5 0
+0.4 0.5 0
+0.4 0.4 0
+0.5 0.4 0
+0.4 0.5 0
+0.5 0.5 0
+0.5 0.4 0
+0.6 0.4 0
+0.5 0.5 0
+0.6 0.5 0
+0.6 0.4 0
+0.7 0.4 0
+0.6 0.5 0
+0.7 0.5 0
+0.7 0.4 0
+0.8 0.4 0
+0.7 0.5 0
+0.8 0.5 0
+0.8 0.4 0
+0.9 0.4 0
+0.8 0.5 0
+0.9 0.5 0
+0.9 0.4 0
+1 0.4 0
+0.9 0.5 0
+1 0.5 0
+0 0.5 0
+0.1 0.5 0
+0 0.6 0
+0.1 0.6 0
+0.1 0.5 0
+0.2 0.5 0
+0.1 0.6 0
+0.2 0.6 0
+0.2 0.5 0
+0.3 0.5 0
+0.2 0.6 0
+0.3 0.6 0
+0.3 0.5 0
+0.4 0.5 0
+0.3 0.6 0
+0.4 0.6 0
+0.4 0.5 0
+0.5 0.5 0
+0.4 0.6 0
+0.5 0.6 0
+0.5 0.5 0
+0.6 0.5 0
+0.5 0.6 0
+0.6 0.6 0
+0.6 0.5 0
+0.7 0.5 0
+0.6 0.6 0
+0.7 0.6 0
+0.7 0.5 0
+0.8 0.5 0
+0.7 0.6 0
+0.8 0.6 0
+0.8 0.5 0
+0.9 0.5 0
+0.8 0.6 0
+0.9 0.6 0
+0.9 0.5 0
+1 0.5 0
+0.9 0.6 0
+1 0.6 0
+0 0.6 0
+0.1 0.6 0
+0 0.7 0
+0.1 0.7 0
+0.1 0.6 0
+0.2 0.6 0
+0.1 0.7 0
+0.2 0.7 0
+0.2 0.6 0
+0.3 0.6 0
+0.2 0.7 0
+0.3 0.7 0
+0.3 0.6 0
+0.4 0.6 0
+0.3 0.7 0
+0.4 0.7 0
+0.4 0.6 0
+0.5 0.6 0
+0.4 0.7 0
+0.5 0.7 0
+0.5 0.6 0
+0.6 0.6 0
+0.5 0.7 0
+0.6 0.7 0
+0.6 0.6 0
+0.7 0.6 0
+0.6 0.7 0
+0.7 0.7 0
+0.7 0.6 0
+0.8 0.6 0
+0.7 0.7 0
+0.8 0.7 0
+0.8 0.6 0
+0.9 0.6 0
+0.8 0.7 0
+0.9 0.7 0
+0.9 0.6 0
+1 0.6 0
+0.9 0.7 0
+1 0.7 0
+0 0.7 0
+0.1 0.7 0
+0 0.8 0
+0.1 0.8 0
+0.1 0.7 0
+0.2 0.7 0
+0.1 0.8 0
+0.2 0.8 0
+0.2 0.7 0
+0.3 0.7 0
+0.2 0.8 0
+0.3 0.8 0
+0.3 0.7 0
+0.4 0.7 0
+0.3 0.8 0
+0.4 0.8 0
+0.4 0.7 0
+0.5 0.7 0
+0.4 0.8 0
+0.5 0.8 0
+0.5 0.7 0
+0.6 0.7 0
+0.5 0.8 0
+0.6 0.8 0
+0.6 0.7 0
+0.7 0.7 0
+0.6 0.8 0
+0.7 0.8 0
+0.7 0.7 0
+0.8 0.7 0
+0.7 0.8 0
+0.8 0.8 0
+0.8 0.7 0
+0.9 0.7 0
+0.8 0.8 0
+0.9 0.8 0
+0.9 0.7 0
+1 0.7 0
+0.9 0.8 0
+1 0.8 0
+0 0.8 0
+0.1 0.8 0
+0 0.9 0
+0.1 0.9 0
+0.1 0.8 0
+0.2 0.8 0
+0.1 0.9 0
+0.2 0.9 0
+0.2 0.8 0
+0.3 0.8 0
+0.2 0.9 0
+0.3 0.9 0
+0.3 0.8 0
+0.4 0.8 0
+0.3 0.9 0
+0.4 0.9 0
+0.4 0.8 0
+0.5 0.8 0
+0.4 0.9 0
+0.5 0.9 0
+0.5 0.8 0
+0.6 0.8 0
+0.5 0.9 0
+0.6 0.9 0
+0.6 0.8 0
+0.7 0.8 0
+0.6 0.9 0
+0.7 0.9 0
+0.7 0.8 0
+0.8 0.8 0
+0.7 0.9 0
+0.8 0.9 0
+0.8 0.8 0
+0.9 0.8 0
+0.8 0.9 0
+0.9 0.9 0
+0.9 0.8 0
+1 0.8 0
+0.9 0.9 0
+1 0.9 0
+0 0.9 0
+0.1 0.9 0
+0 1 0
+0.1 1 0
+0.1 0.9 0
+0.2 0.9 0
+0.1 1 0
+0.2 1 0
+0.2 0.9 0
+0.3 0.9 0
+0.2 1 0
+0.3 1 0
+0.3 0.9 0
+0.4 0.9 0
+0.3 1 0
+0.4 1 0
+0.4 0.9 0
+0.5 0.9 0
+0.4 1 0
+0.5 1 0
+0.5 0.9 0
+0.6 0.9 0
+0.5 1 0
+0.6 1 0
+0.6 0.9 0
+0.7 0.9 0
+0.6 1 0
+0.7 1 0
+0.7 0.9 0
+0.8 0.9 0
+0.7 1 0
+0.8 1 0
+0.8 0.9 0
+0.9 0.9 0
+0.8 1 0
+0.9 1 0
+0.9 0.9 0
+1 0.9 0
+0.9 1 0
+1 1 0
+</DataArray>
+</Points>
+<Cells>
+<DataArray type="Int32" Name="connectivity" format="ascii">
+0 1 3 2 
+4 5 7 6 
+8 9 11 10 
+12 13 15 14 
+16 17 19 18 
+20 21 23 22 
+24 25 27 26 
+28 29 31 30 
+32 33 35 34 
+36 37 39 38 
+40 41 43 42 
+44 45 47 46 
+48 49 51 50 
+52 53 55 54 
+56 57 59 58 
+60 61 63 62 
+64 65 67 66 
+68 69 71 70 
+72 73 75 74 
+76 77 79 78 
+80 81 83 82 
+84 85 87 86 
+88 89 91 90 
+92 93 95 94 
+96 97 99 98 
+100 101 103 102 
+104 105 107 106 
+108 109 111 110 
+112 113 115 114 
+116 117 119 118 
+120 121 123 122 
+124 125 127 126 
+128 129 131 130 
+132 133 135 134 
+136 137 139 138 
+140 141 143 142 
+144 145 147 146 
+148 149 151 150 
+152 153 155 154 
+156 157 159 158 
+160 161 163 162 
+164 165 167 166 
+168 169 171 170 
+172 173 175 174 
+176 177 179 178 
+180 181 183 182 
+184 185 187 186 
+188 189 191 190 
+192 193 195 194 
+196 197 199 198 
+200 201 203 202 
+204 205 207 206 
+208 209 211 210 
+212 213 215 214 
+216 217 219 218 
+220 221 223 222 
+224 225 227 226 
+228 229 231 230 
+232 233 235 234 
+236 237 239 238 
+240 241 243 242 
+244 245 247 246 
+248 249 251 250 
+252 253 255 254 
+256 257 259 258 
+260 261 263 262 
+264 265 267 266 
+268 269 271 270 
+272 273 275 274 
+276 277 279 278 
+280 281 283 282 
+284 285 287 286 
+288 289 291 290 
+292 293 295 294 
+296 297 299 298 
+300 301 303 302 
+304 305 307 306 
+308 309 311 310 
+312 313 315 314 
+316 317 319 318 
+320 321 323 322 
+324 325 327 326 
+328 329 331 330 
+332 333 335 334 
+336 337 339 338 
+340 341 343 342 
+344 345 347 346 
+348 349 351 350 
+352 353 355 354 
+356 357 359 358 
+360 361 363 362 
+364 365 367 366 
+368 369 371 370 
+372 373 375 374 
+376 377 379 378 
+380 381 383 382 
+384 385 387 386 
+388 389 391 390 
+392 393 395 394 
+396 397 399 398 
+</DataArray>
+<DataArray type="Int32" Name="offsets" format="ascii">
+4
+8
+12
+16
+20
+24
+28
+32
+36
+40
+44
+48
+52
+56
+60
+64
+68
+72
+76
+80
+84
+88
+92
+96
+100
+104
+108
+112
+116
+120
+124
+128
+132
+136
+140
+144
+148
+152
+156
+160
+164
+168
+172
+176
+180
+184
+188
+192
+196
+200
+204
+208
+212
+216
+220
+224
+228
+232
+236
+240
+244
+248
+252
+256
+260
+264
+268
+272
+276
+280
+284
+288
+292
+296
+300
+304
+308
+312
+316
+320
+324
+328
+332
+336
+340
+344
+348
+352
+356
+360
+364
+368
+372
+376
+380
+384
+388
+392
+396
+400
+</DataArray>
+<DataArray type="UInt8" Name="types" format="ascii">
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+9
+</DataArray>
+</Cells>
+<CellData Scalars="attribute">
+<DataArray type="Int32" Name="attribute" format="ascii">
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+</DataArray>
+</CellData>
+<PointData >
+<DataArray type="Float64" Name="concentration" NumberOfComponents="1"  format="ascii" >
+1
+1
+1.34164079
+1.34164079
+1
+1
+1.34164079
+1.34164079
+1
+1
+1.34164079
+1.34164079
+1
+1
+1.34164079
+1.34164079
+1
+1
+1.34164079
+1.34164079
+1
+1
+1.34164079
+1.34164079
+1
+1
+1.34164079
+1.34164079
+1
+1
+1.34164079
+1.34164079
+1
+1
+1.34164079
+1.34164079
+1
+1
+1.34164079
+1.34164079
+1.34164079
+1.34164079
+1.61245155
+1.61245155
+1.34164079
+1.34164079
+1.61245155
+1.61245155
+1.34164079
+1.34164079
+1.61245155
+1.61245155
+1.34164079
+1.34164079
+1.61245155
+1.61245155
+1.34164079
+1.34164079
+1.61245155
+1.61245155
+1.34164079
+1.34164079
+1.61245155
+1.61245155
+1.34164079
+1.34164079
+1.61245155
+1.61245155
+1.34164079
+1.34164079
+1.61245155
+1.61245155
+1.34164079
+1.34164079
+1.61245155
+1.61245155
+1.34164079
+1.34164079
+1.61245155
+1.61245155
+1.61245155
+1.61245155
+1.84390889
+1.84390889
+1.61245155
+1.61245155
+1.84390889
+1.84390889
+1.61245155
+1.61245155
+1.84390889
+1.84390889
+1.61245155
+1.61245155
+1.84390889
+1.84390889
+1.61245155
+1.61245155
+1.84390889
+1.84390889
+1.61245155
+1.61245155
+1.84390889
+1.84390889
+1.61245155
+1.61245155
+1.84390889
+1.84390889
+1.61245155
+1.61245155
+1.84390889
+1.84390889
+1.61245155
+1.61245155
+1.84390889
+1.84390889
+1.61245155
+1.61245155
+1.84390889
+1.84390889
+1.84390889
+1.84390889
+2.04939015
+2.04939015
+1.84390889
+1.84390889
+2.04939015
+2.04939015
+1.84390889
+1.84390889
+2.04939015
+2.04939015
+1.84390889
+1.84390889
+2.04939015
+2.04939015
+1.84390889
+1.84390889
+2.04939015
+2.04939015
+1.84390889
+1.84390889
+2.04939015
+2.04939015
+1.84390889
+1.84390889
+2.04939015
+2.04939015
+1.84390889
+1.84390889
+2.04939015
+2.04939015
+1.84390889
+1.84390889
+2.04939015
+2.04939015
+1.84390889
+1.84390889
+2.04939015
+2.04939015
+2.04939015
+2.04939015
+2.23606798
+2.23606798
+2.04939015
+2.04939015
+2.23606798
+2.23606798
+2.04939015
+2.04939015
+2.23606798
+2.23606798
+2.04939015
+2.04939015
+2.23606798
+2.23606798
+2.04939015
+2.04939015
+2.23606798
+2.23606798
+2.04939015
+2.04939015
+2.23606798
+2.23606798
+2.04939015
+2.04939015
+2.23606798
+2.23606798
+2.04939015
+2.04939015
+2.23606798
+2.23606798
+2.04939015
+2.04939015
+2.23606798
+2.23606798
+2.04939015
+2.04939015
+2.23606798
+2.23606798
+2.23606798
+2.23606798
+2.40831892
+2.40831892
+2.23606798
+2.23606798
+2.40831892
+2.40831892
+2.23606798
+2.23606798
+2.40831892
+2.40831892
+2.23606798
+2.23606798
+2.40831892
+2.40831892
+2.23606798
+2.23606798
+2.40831892
+2.40831892
+2.23606798
+2.23606798
+2.40831892
+2.40831892
+2.23606798
+2.23606798
+2.40831892
+2.40831892
+2.23606798
+2.23606798
+2.40831892
+2.40831892
+2.23606798
+2.23606798
+2.40831892
+2.40831892
+2.23606798
+2.23606798
+2.40831892
+2.40831892
+2.40831892
+2.40831892
+2.56904652
+2.56904652
+2.40831892
+2.40831892
+2.56904652
+2.56904652
+2.40831892
+2.40831892
+2.56904652
+2.56904652
+2.40831892
+2.40831892
+2.56904652
+2.56904652
+2.40831892
+2.40831892
+2.56904652
+2.56904652
+2.40831892
+2.40831892
+2.56904652
+2.56904652
+2.40831892
+2.40831892
+2.56904652
+2.56904652
+2.40831892
+2.40831892
+2.56904652
+2.56904652
+2.40831892
+2.40831892
+2.56904652
+2.56904652
+2.40831892
+2.40831892
+2.56904652
+2.56904652
+2.56904652
+2.56904652
+2.7202941
+2.7202941
+2.56904652
+2.56904652
+2.7202941
+2.7202941
+2.56904652
+2.56904652
+2.7202941
+2.7202941
+2.56904652
+2.56904652
+2.7202941
+2.7202941
+2.56904652
+2.56904652
+2.7202941
+2.7202941
+2.56904652
+2.56904652
+2.7202941
+2.7202941
+2.56904652
+2.56904652
+2.7202941
+2.7202941
+2.56904652
+2.56904652
+2.7202941
+2.7202941
+2.56904652
+2.56904652
+2.7202941
+2.7202941
+2.56904652
+2.56904652
+2.7202941
+2.7202941
+2.7202941
+2.7202941
+2.86356421
+2.86356421
+2.7202941
+2.7202941
+2.86356421
+2.86356421
+2.7202941
+2.7202941
+2.86356421
+2.86356421
+2.7202941
+2.7202941
+2.86356421
+2.86356421
+2.7202941
+2.7202941
+2.86356421
+2.86356421
+2.7202941
+2.7202941
+2.86356421
+2.86356421
+2.7202941
+2.7202941
+2.86356421
+2.86356421
+2.7202941
+2.7202941
+2.86356421
+2.86356421
+2.7202941
+2.7202941
+2.86356421
+2.86356421
+2.7202941
+2.7202941
+2.86356421
+2.86356421
+2.86356421
+2.86356421
+3
+3
+2.86356421
+2.86356421
+3
+3
+2.86356421
+2.86356421
+3
+3
+2.86356421
+2.86356421
+3
+3
+2.86356421
+2.86356421
+3
+3
+2.86356421
+2.86356421
+3
+3
+2.86356421
+2.86356421
+3
+3
+2.86356421
+2.86356421
+3
+3
+2.86356421
+2.86356421
+3
+3
+2.86356421
+2.86356421
+3
+3
+</DataArray>
+</PointData>
+</Piece>
+</UnstructuredGrid>
+</VTKFile>
diff --git a/test/tests/mfem/kernels/gold/OutputData/NLDiffusion/Run0/Run0.pvd b/test/tests/mfem/kernels/gold/OutputData/NLDiffusion/Run0/Run0.pvd
new file mode 100644
index 000000000000..72ed1b22e1dd
--- /dev/null
+++ b/test/tests/mfem/kernels/gold/OutputData/NLDiffusion/Run0/Run0.pvd
@@ -0,0 +1,7 @@
+<?xml version="1.0"?>
+<VTKFile type="Collection" version="2.2" byte_order="LittleEndian">
+<Collection>
+<DataSet timestep="0" group="" part="0" file="Cycle000000/data.pvtu" name="mesh"/>
+<DataSet timestep="1" group="" part="0" file="Cycle000001/data.pvtu" name="mesh"/>
+</Collection>
+</VTKFile>
diff --git a/test/tests/mfem/kernels/nldiffusion.i b/test/tests/mfem/kernels/nldiffusion.i
new file mode 100644
index 000000000000..36a1212a4269
--- /dev/null
+++ b/test/tests/mfem/kernels/nldiffusion.i
@@ -0,0 +1,111 @@
+[Mesh]
+  type = MFEMMesh
+  file = ../mesh/square.e
+  dim = 2
+[]
+
+[Problem]
+  type = MFEMProblem
+[]
+
+[FESpaces]
+  [H1FESpace]
+    type = MFEMScalarFESpace
+    fec_type = H1
+    fec_order = FIRST
+  []
+[]
+
+[Variables]
+  [concentration]
+    type = MFEMVariable
+    fespace = H1FESpace
+  []
+[]
+
+
+[ICs]
+  [diffused_ic]
+    type = MFEMScalarIC
+    coefficient = initial
+    variable = concentration
+  []
+[]
+
+[Functions]
+  [one]
+    type = ParsedFunction
+    expression = 1.0
+  []
+  [initial]
+    type = ParsedFunction
+    expression = 2*y+1
+  []
+  [diffCoeff]
+    type = MFEMParsedFunction
+    expression = u
+    symbol_names = 'u'
+    symbol_values = 'concentration'
+  []
+[]
+
+[BCs]
+  [top]
+    type = MFEMScalarDirichletBC
+    variable = concentration
+    boundary = 'top'
+    coefficient = 3.0
+  []
+  [bottom]
+    type = MFEMScalarDirichletBC
+    variable = concentration
+    boundary = 'bottom'
+    coefficient = 1.0
+  []
+
+[]
+
+[Kernels]
+ active = 'nl'
+  [nl]
+    type = MFEMNLDiffusionKernel
+    variable = concentration
+    coefficient = diffCoeff
+  []
+[]
+
+[Preconditioner]
+  [boomeramg]
+    type = MFEMHypreBoomerAMG
+    print_level = 0
+  []
+  [jacobi]
+    type = MFEMOperatorJacobiSmoother
+  []
+[]
+
+[Solver]
+  type = MFEMHypreGMRES
+  preconditioner = boomeramg
+  print_level = 1
+  l_tol = 1e-16
+  l_max_its = 1000
+[]
+
+[Executioner]
+  type = MFEMSteady
+  device = cpu
+  nl_max_its = 100
+  nl_abs_tol = 1.0e-10
+  nl_rel_tol = 1.0e-9
+  print_level = 1
+[]
+
+[Outputs]
+  active = ParaViewDataCollection
+  [ParaViewDataCollection]
+    type = MFEMParaViewDataCollection
+    file_base = OutputData/NLDiffusion
+    vtk_format = ASCII
+  []
+[]
diff --git a/test/tests/mfem/kernels/tests b/test/tests/mfem/kernels/tests
index 07ce96e449d6..46008c91749c 100644
--- a/test/tests/mfem/kernels/tests
+++ b/test/tests/mfem/kernels/tests
@@ -188,4 +188,15 @@
     compute_devices = 'cpu cuda'
     recover = false
   []
+  [MFEMNonLinearDiffusion]
+    type = XMLDiff
+    input = nldiffusion.i
+    xmldiff = 'OutputData/NLDiffusion/Run0/Run0.pvd
+                OutputData/NLDiffusion/Run0/Cycle000001/proc000000.vtu'
+    requirement = 'The system shall have the ability to solve a non-linear diffusion problem using MFEM.'
+    capabilities = 'mfem'
+    compute_devices = 'cpu'
+    max_parallel = 1 # schemadiff with multiple ranks
+    recover = false
+  []  
 []