diff --git a/framework/include/mfem/bcs/MFEMBoundaryCondition.h b/framework/include/mfem/bcs/MFEMBoundaryCondition.h
index 75af5a8c..c25d3473 100644
--- a/framework/include/mfem/bcs/MFEMBoundaryCondition.h
+++ b/framework/include/mfem/bcs/MFEMBoundaryCondition.h
@@ -1,9 +1,9 @@
 #pragma once
 
 #include "GeneralUserObject.h"
-#include "boundary_conditions.hpp"
-#include "gridfunctions.hpp"
-#include "coefficients.hpp"
+#include "boundary_conditions.h"
+#include "gridfunctions.h"
+#include "coefficients.h"
 #include "Function.h"
 
 class MFEMBoundaryCondition : public GeneralUserObject
@@ -14,7 +14,7 @@ class MFEMBoundaryCondition : public GeneralUserObject
   MFEMBoundaryCondition(const InputParameters & parameters);
   ~MFEMBoundaryCondition() override {}
 
-  inline virtual std::shared_ptr<hephaestus::BoundaryCondition> getBC() const
+  inline virtual std::shared_ptr<platypus::BoundaryCondition> getBC() const
   {
     return _boundary_condition;
   }
@@ -27,5 +27,5 @@ class MFEMBoundaryCondition : public GeneralUserObject
   std::vector<BoundaryName> _boundary_names;
   mfem::Array<int> bdr_attr;
 
-  std::shared_ptr<hephaestus::BoundaryCondition> _boundary_condition{nullptr};
+  std::shared_ptr<platypus::BoundaryCondition> _boundary_condition{nullptr};
 };
diff --git a/framework/include/mfem/bcs/MFEMScalarDirichletBC.h b/framework/include/mfem/bcs/MFEMScalarDirichletBC.h
index 3455249b..9ff404cc 100644
--- a/framework/include/mfem/bcs/MFEMScalarDirichletBC.h
+++ b/framework/include/mfem/bcs/MFEMScalarDirichletBC.h
@@ -2,7 +2,7 @@
 
 #include "MFEMBoundaryCondition.h"
 #include "MFEMFunctionCoefficient.h"
-#include "boundary_conditions.hpp"
+#include "boundary_conditions.h"
 
 class MFEMScalarDirichletBC : public MFEMBoundaryCondition
 {
diff --git a/framework/include/mfem/bcs/MFEMVectorDirichletBC.h b/framework/include/mfem/bcs/MFEMVectorDirichletBC.h
index cfab0b1d..32232e78 100644
--- a/framework/include/mfem/bcs/MFEMVectorDirichletBC.h
+++ b/framework/include/mfem/bcs/MFEMVectorDirichletBC.h
@@ -2,7 +2,7 @@
 
 #include "MFEMBoundaryCondition.h"
 #include "MFEMVectorFunctionCoefficient.h"
-#include "boundary_conditions.hpp"
+#include "boundary_conditions.h"
 
 class MFEMVectorDirichletBC : public MFEMBoundaryCondition
 {
diff --git a/framework/include/mfem/bcs/MFEMVectorNormalIntegratedBC.h b/framework/include/mfem/bcs/MFEMVectorNormalIntegratedBC.h
index 4140d56e..08970b56 100644
--- a/framework/include/mfem/bcs/MFEMVectorNormalIntegratedBC.h
+++ b/framework/include/mfem/bcs/MFEMVectorNormalIntegratedBC.h
@@ -2,7 +2,7 @@
 
 #include "MFEMBoundaryCondition.h"
 #include "MFEMVectorFunctionCoefficient.h"
-#include "boundary_conditions.hpp"
+#include "boundary_conditions.h"
 
 class MFEMVectorNormalIntegratedBC : public MFEMBoundaryCondition
 {
diff --git a/framework/include/mfem/bcs/boundary_conditions.h b/framework/include/mfem/bcs/boundary_conditions.h
new file mode 100644
index 00000000..18f88174
--- /dev/null
+++ b/framework/include/mfem/bcs/boundary_conditions.h
@@ -0,0 +1,36 @@
+#pragma once
+#include "essential_bcs.h"
+#include "integrated_bcs.h"
+#include "named_fields_map.h"
+#include "robin_bcs.h"
+
+namespace platypus
+{
+
+class BCMap : public platypus::NamedFieldsMap<platypus::BoundaryCondition>
+{
+public:
+  mfem::Array<int> GetEssentialBdrMarkers(const std::string & name_, mfem::Mesh * mesh_);
+
+  void ApplyEssentialBCs(const std::string & name_,
+                         mfem::Array<int> & ess_tdof_list,
+                         mfem::GridFunction & gridfunc,
+                         mfem::Mesh * mesh_);
+
+  void ApplyEssentialBCs(const std::string & name_,
+                         mfem::Array<int> & ess_tdof_list,
+                         mfem::ParComplexGridFunction & gridfunc,
+                         mfem::Mesh * mesh_);
+
+  void ApplyIntegratedBCs(const std::string & name_, mfem::LinearForm & lf, mfem::Mesh * mesh_);
+
+  void ApplyIntegratedBCs(const std::string & name_,
+                          mfem::ParComplexLinearForm & clf,
+                          mfem::Mesh * mesh_);
+
+  void ApplyIntegratedBCs(const std::string & name_,
+                          mfem::ParSesquilinearForm & clf,
+                          mfem::Mesh * mesh_);
+};
+
+} // namespace platypus
diff --git a/framework/include/mfem/equation_systems/equation_system.h b/framework/include/mfem/equation_systems/equation_system.h
new file mode 100644
index 00000000..0f380ead
--- /dev/null
+++ b/framework/include/mfem/equation_systems/equation_system.h
@@ -0,0 +1,141 @@
+#pragma once
+#include "../common/pfem_extras.hpp"
+#include "inputs.h"
+#include "kernel_base.h"
+#include "named_fields_map.h"
+
+namespace platypus
+{
+
+/*
+Class to store weak form components (bilinear and linear forms, and optionally
+mixed and nonlinear forms) and build methods
+*/
+class EquationSystem : public mfem::Operator
+{
+public:
+  using ParBilinearFormKernel = platypus::Kernel<mfem::ParBilinearForm>;
+  using ParLinearFormKernel = platypus::Kernel<mfem::ParLinearForm>;
+  using ParNonlinearFormKernel = platypus::Kernel<mfem::ParNonlinearForm>;
+  using ParMixedBilinearFormKernel = platypus::Kernel<mfem::ParMixedBilinearForm>;
+
+  EquationSystem() = default;
+  ~EquationSystem() override;
+
+  // Test variables are associated with LinearForms,
+  // whereas trial variables are associated with gridfunctions.
+
+  // Names of all variables corresponding to gridfunctions. This may differ
+  // from test_var_names when time derivatives are present.
+  std::vector<std::string> _trial_var_names;
+  // Names of all test variables corresponding to linear forms in this equation
+  // system
+  std::vector<std::string> _test_var_names;
+  std::vector<mfem::ParFiniteElementSpace *> _test_pfespaces;
+
+  // Components of weak form. // Named according to test variable
+  platypus::NamedFieldsMap<mfem::ParBilinearForm> _blfs;
+  platypus::NamedFieldsMap<mfem::ParLinearForm> _lfs;
+  platypus::NamedFieldsMap<mfem::ParNonlinearForm> _nlfs;
+  platypus::NamedFieldsMap<platypus::NamedFieldsMap<mfem::ParMixedBilinearForm>>
+      _mblfs; // named according to trial variable
+
+  // add test variable to EquationSystem;
+  virtual void AddTestVariableNameIfMissing(const std::string & test_var_name);
+  virtual void AddTrialVariableNameIfMissing(const std::string & trial_var_name);
+
+  // Add kernels.
+  void AddKernel(const std::string & test_var_name,
+                 std::shared_ptr<ParBilinearFormKernel> blf_kernel);
+
+  void AddKernel(const std::string & test_var_name, std::shared_ptr<ParLinearFormKernel> lf_kernel);
+
+  void AddKernel(const std::string & test_var_name,
+                 std::shared_ptr<ParNonlinearFormKernel> nlf_kernel);
+
+  void AddKernel(const std::string & trial_var_name,
+                 const std::string & test_var_name,
+                 std::shared_ptr<ParMixedBilinearFormKernel> mblf_kernel);
+
+  virtual void ApplyBoundaryConditions(platypus::BCMap & bc_map);
+
+  // override to add kernels
+  virtual void AddKernels() {}
+
+  // Build forms
+  virtual void Init(platypus::GridFunctions & gridfunctions,
+                    const platypus::FESpaces & fespaces,
+                    platypus::BCMap & bc_map,
+                    platypus::Coefficients & coefficients);
+  virtual void BuildLinearForms(platypus::BCMap & bc_map);
+  virtual void BuildBilinearForms();
+  virtual void BuildMixedBilinearForms();
+  virtual void BuildEquationSystem(platypus::BCMap & bc_map);
+
+  // Form linear system, with essential boundary conditions accounted for
+  virtual void FormLinearSystem(mfem::OperatorHandle & op,
+                                mfem::BlockVector & trueX,
+                                mfem::BlockVector & trueRHS);
+
+  // Build linear system, with essential boundary conditions accounted for
+  virtual void BuildJacobian(mfem::BlockVector & trueX, mfem::BlockVector & trueRHS);
+
+  /// Compute residual y = Mu
+  void Mult(const mfem::Vector & u, mfem::Vector & residual) const override;
+
+  /// Compute J = M + grad_H(u)
+  mfem::Operator & GetGradient(const mfem::Vector & u) const override;
+
+  // Update variable from solution vector after solve
+  virtual void RecoverFEMSolution(mfem::BlockVector & trueX,
+                                  platypus::GridFunctions & gridfunctions);
+
+  std::vector<mfem::Array<int>> _ess_tdof_lists;
+
+protected:
+  bool VectorContainsName(const std::vector<std::string> & the_vector,
+                          const std::string & name) const;
+
+  // gridfunctions for setting Dirichlet BCs
+  std::vector<std::unique_ptr<mfem::ParGridFunction>> _xs;
+
+  mfem::Array2D<mfem::HypreParMatrix *> _h_blocks;
+
+  // Arrays to store kernels to act on each component of weak form. Named
+  // according to test variable
+  platypus::NamedFieldsMap<std::vector<std::shared_ptr<ParBilinearFormKernel>>> _blf_kernels_map;
+
+  platypus::NamedFieldsMap<std::vector<std::shared_ptr<ParLinearFormKernel>>> _lf_kernels_map;
+
+  platypus::NamedFieldsMap<std::vector<std::shared_ptr<ParNonlinearFormKernel>>> _nlf_kernels_map;
+
+  platypus::NamedFieldsMap<
+      platypus::NamedFieldsMap<std::vector<std::shared_ptr<ParMixedBilinearFormKernel>>>>
+      _mblf_kernels_map_map;
+
+  mutable mfem::OperatorHandle _jacobian;
+};
+
+/*
+Class to store weak form components for time dependent PDEs
+*/
+class TimeDependentEquationSystem : public EquationSystem
+{
+public:
+  TimeDependentEquationSystem();
+  ~TimeDependentEquationSystem() override = default;
+
+  static std::string GetTimeDerivativeName(std::string name)
+  {
+    return std::string("d") + name + std::string("_dt");
+  }
+
+  void AddTrialVariableNameIfMissing(const std::string & trial_var_name) override;
+
+  virtual void SetTimeStep(double dt);
+  virtual void UpdateEquationSystem(platypus::BCMap & bc_map);
+  mfem::ConstantCoefficient _dt_coef; // Coefficient for timestep scaling
+  std::vector<std::string> _trial_var_time_derivative_names;
+};
+
+} // namespace platypus
diff --git a/framework/include/mfem/fespaces/MFEMFESpace.h b/framework/include/mfem/fespaces/MFEMFESpace.h
index 0e01ae91..67a5e9b0 100644
--- a/framework/include/mfem/fespaces/MFEMFESpace.h
+++ b/framework/include/mfem/fespaces/MFEMFESpace.h
@@ -1,8 +1,5 @@
 #pragma once
-
 #include "GeneralUserObject.h"
-#include "inputs.hpp"
-#include "gridfunctions.hpp"
 
 class MFEMFESpace : public GeneralUserObject
 {
diff --git a/framework/include/mfem/kernels/MFEMDiffusionKernel.h b/framework/include/mfem/kernels/MFEMDiffusionKernel.h
index 4d75fd62..b8c3cf6d 100644
--- a/framework/include/mfem/kernels/MFEMDiffusionKernel.h
+++ b/framework/include/mfem/kernels/MFEMDiffusionKernel.h
@@ -1,6 +1,6 @@
 #pragma once
 #include "MFEMBilinearFormKernel.h"
-#include "kernels.hpp"
+#include "kernels.h"
 
 class MFEMDiffusionKernel : public MFEMBilinearFormKernel
 {
@@ -14,12 +14,9 @@ class MFEMDiffusionKernel : public MFEMBilinearFormKernel
   virtual void initialize() override {}
   virtual void finalize() override {}
 
-  std::shared_ptr<hephaestus::Kernel<mfem::ParBilinearForm>> getKernel() override
-  {
-    return _kernel;
-  }
+  std::shared_ptr<platypus::Kernel<mfem::ParBilinearForm>> getKernel() override { return _kernel; }
 
 protected:
-  hephaestus::InputParameters _kernel_params;
-  std::shared_ptr<hephaestus::DiffusionKernel> _kernel{nullptr};
+  platypus::InputParameters _kernel_params;
+  std::shared_ptr<platypus::DiffusionKernel> _kernel{nullptr};
 };
diff --git a/framework/include/mfem/materials/MFEMMaterial.h b/framework/include/mfem/materials/MFEMMaterial.h
index efadda0c..d4a3e5c5 100644
--- a/framework/include/mfem/materials/MFEMMaterial.h
+++ b/framework/include/mfem/materials/MFEMMaterial.h
@@ -2,7 +2,7 @@
 
 #include "GeneralUserObject.h"
 #include "MFEMCoefficient.h"
-#include "coefficients.hpp"
+#include "coefficients.h"
 
 class MFEMMaterial : public GeneralUserObject
 {
@@ -16,7 +16,7 @@ class MFEMMaterial : public GeneralUserObject
   virtual void initialize() override {}
   virtual void finalize() override {}
 
-  virtual void storeCoefficients(hephaestus::Subdomain & subdomain) {}
+  virtual void storeCoefficients(platypus::Subdomain & subdomain) {}
 
   std::vector<SubdomainName> blocks;
 };
diff --git a/framework/include/mfem/problem/MFEMProblem.h b/framework/include/mfem/problem/MFEMProblem.h
index 1be63adc..4481413d 100644
--- a/framework/include/mfem/problem/MFEMProblem.h
+++ b/framework/include/mfem/problem/MFEMProblem.h
@@ -15,7 +15,6 @@
 #include "MFEMBilinearFormKernel.h"
 #include "MFEMLinearFormKernel.h"
 #include "MFEMFormulation.h"
-#include "MFEMAuxSolver.h"
 #include "MFEMDataCollection.h"
 #include "MFEMFESpace.h"
 #include "Function.h"
@@ -23,7 +22,7 @@
 #include "SystemBase.h"
 #include "Transient.h"
 #include "Steady.h"
-#include "hephaestus.hpp"
+#include "platypus.h"
 #include "libmesh/string_to_enum.h"
 #include "libmesh/point.h"
 
@@ -109,9 +108,7 @@ class MFEMProblem : public ExternalProblem
                  const std::string & name,
                  InputParameters & parameters) override;
   /**
-   * Override of ExternalProblem::addAuxKernel. Uses ExternalProblem::addAuxKernel to create a
-   * GeneralUserObject representing the auxkernel in MOOSE, and creates corresponding MFEM auxsolver
-   * to be used in the MFEM solve.
+   * Override of ExternalProblem::addAuxKernel.
    */
   void addAuxKernel(const std::string & kernel_name,
                     const std::string & name,
@@ -151,9 +148,9 @@ class MFEMProblem : public ExternalProblem
    * builders.
    */
   template <class T>
-  void addKernel(std::string var_name, std::shared_ptr<hephaestus::Kernel<T>> kernel)
+  void addKernel(std::string var_name, std::shared_ptr<platypus::Kernel<T>> kernel)
   {
-    using namespace hephaestus;
+    using namespace platypus;
 
     EquationSystemProblemBuilderInterface * eqn_system_problem_builder{nullptr};
 
@@ -173,13 +170,13 @@ class MFEMProblem : public ExternalProblem
   std::string _formulation_name;
   int _order;
 
-  hephaestus::Coefficients _coefficients;
-  hephaestus::InputParameters _solver_options;
-  hephaestus::Outputs _outputs;
-  hephaestus::InputParameters _exec_params;
+  platypus::Coefficients _coefficients;
+  platypus::InputParameters _solver_options;
+  platypus::Outputs _outputs;
+  platypus::InputParameters _exec_params;
 
-  std::shared_ptr<hephaestus::ProblemBuilder> mfem_problem_builder{nullptr};
+  std::shared_ptr<platypus::ProblemBuilder> mfem_problem_builder{nullptr};
 
-  std::unique_ptr<hephaestus::Problem> mfem_problem{nullptr};
-  std::unique_ptr<hephaestus::Executioner> executioner{nullptr};
+  std::unique_ptr<platypus::Problem> mfem_problem{nullptr};
+  std::unique_ptr<platypus::Executioner> executioner{nullptr};
 };
diff --git a/framework/include/mfem/problem_operators/equation_system_problem_operator.h b/framework/include/mfem/problem_operators/equation_system_problem_operator.h
new file mode 100644
index 00000000..9b507a18
--- /dev/null
+++ b/framework/include/mfem/problem_operators/equation_system_problem_operator.h
@@ -0,0 +1,39 @@
+#pragma once
+#include "problem_operator.h"
+#include "problem_operator_interface.h"
+#include "equation_system_interface.h"
+
+namespace platypus
+{
+/// Steady-state problem operator with an equation system.
+class EquationSystemProblemOperator : public ProblemOperator, public EquationSystemInterface
+{
+public:
+  EquationSystemProblemOperator(platypus::Problem &) = delete;
+
+  EquationSystemProblemOperator(platypus::Problem & problem,
+                                std::unique_ptr<platypus::EquationSystem> equation_system)
+    : ProblemOperator(problem), _equation_system(std::move(equation_system))
+  {
+  }
+
+  void SetGridFunctions() override;
+  void Init(mfem::Vector & X) override;
+
+  ~EquationSystemProblemOperator() override = default;
+
+  [[nodiscard]] platypus::EquationSystem * GetEquationSystem() const override
+  {
+    if (!_equation_system)
+    {
+      MFEM_ABORT("No equation system has been added to ProblemOperator.");
+    }
+
+    return _equation_system.get();
+  }
+
+private:
+  std::unique_ptr<platypus::EquationSystem> _equation_system{nullptr};
+};
+
+} // namespace platypus
\ No newline at end of file
diff --git a/framework/include/mfem/problem_operators/problem_operator.h b/framework/include/mfem/problem_operators/problem_operator.h
new file mode 100644
index 00000000..f07c3947
--- /dev/null
+++ b/framework/include/mfem/problem_operators/problem_operator.h
@@ -0,0 +1,22 @@
+#pragma once
+#include "../common/pfem_extras.hpp"
+#include "hephaestus_solvers.h"
+#include "problem_builder_base.h"
+#include "problem_operator_interface.h"
+
+namespace platypus
+{
+/// Steady-state problem operator with no equation system.
+class ProblemOperator : public mfem::Operator, public ProblemOperatorInterface
+{
+public:
+  ProblemOperator(platypus::Problem & problem) : ProblemOperatorInterface(problem) {}
+  ~ProblemOperator() override = default;
+
+  void SetGridFunctions() override;
+
+  virtual void Solve(mfem::Vector & X) {}
+  void Mult(const mfem::Vector & x, mfem::Vector & y) const override {}
+};
+
+} // namespace platypus
\ No newline at end of file
diff --git a/framework/include/mfem/problem_operators/problem_operator_interface.h b/framework/include/mfem/problem_operators/problem_operator_interface.h
new file mode 100644
index 00000000..efd3ab88
--- /dev/null
+++ b/framework/include/mfem/problem_operators/problem_operator_interface.h
@@ -0,0 +1,30 @@
+#pragma once
+#include "problem_builder_base.h"
+
+namespace platypus
+{
+/// Interface inherited by ProblemOperator and TimeDomainProblemOperator. Removes duplicated code in both classes.
+class ProblemOperatorInterface
+{
+public:
+  ProblemOperatorInterface(platypus::Problem & problem) : _problem(problem) {}
+  virtual ~ProblemOperatorInterface() = default;
+
+  virtual void SetGridFunctions();
+  virtual void Init(mfem::Vector & X);
+
+  mfem::Array<int> _true_offsets, _block_true_offsets;
+
+  mfem::BlockVector _true_x, _true_rhs;
+  mfem::OperatorHandle _equation_system_operator;
+
+protected:
+  // Reference to the current problem.
+  platypus::Problem & _problem;
+
+  // Vector of names of state gridfunctions used in formulation, ordered by appearance in block
+  // vector during solve.
+  std::vector<std::string> _trial_var_names;
+  std::vector<mfem::ParGridFunction *> _trial_variables;
+};
+}
\ No newline at end of file
diff --git a/framework/include/mfem/problem_operators/time_domain_equation_system_problem_operator.h b/framework/include/mfem/problem_operators/time_domain_equation_system_problem_operator.h
new file mode 100644
index 00000000..970aa421
--- /dev/null
+++ b/framework/include/mfem/problem_operators/time_domain_equation_system_problem_operator.h
@@ -0,0 +1,46 @@
+#pragma once
+#include "../common/pfem_extras.hpp"
+#include "time_domain_problem_operator.h"
+#include "problem_operator_interface.h"
+#include "equation_system_interface.h"
+
+namespace platypus
+{
+
+/// Problem operator for time-dependent problems with an equation system.
+class TimeDomainEquationSystemProblemOperator : public TimeDomainProblemOperator,
+                                                public EquationSystemInterface
+{
+public:
+  TimeDomainEquationSystemProblemOperator(platypus::Problem &) = delete;
+  TimeDomainEquationSystemProblemOperator(
+      platypus::Problem & problem,
+      std::unique_ptr<platypus::TimeDependentEquationSystem> equation_system)
+    : TimeDomainProblemOperator(problem), _equation_system{std::move(equation_system)}
+  {
+  }
+
+  void SetGridFunctions() override;
+  void Init(mfem::Vector & X) override;
+
+  void ImplicitSolve(const double dt, const mfem::Vector & X, mfem::Vector & dX_dt) override;
+
+  [[nodiscard]] platypus::TimeDependentEquationSystem * GetEquationSystem() const override
+  {
+    if (!_equation_system)
+    {
+      MFEM_ABORT("No equation system has been added.");
+    }
+
+    return _equation_system.get();
+  }
+
+protected:
+  void BuildEquationSystemOperator(double dt);
+
+private:
+  std::vector<mfem::ParGridFunction *> _trial_variable_time_derivatives;
+  std::unique_ptr<platypus::TimeDependentEquationSystem> _equation_system{nullptr};
+};
+
+} // namespace platypus
\ No newline at end of file
diff --git a/framework/include/mfem/problem_operators/time_domain_problem_operator.h b/framework/include/mfem/problem_operators/time_domain_problem_operator.h
new file mode 100644
index 00000000..73160f44
--- /dev/null
+++ b/framework/include/mfem/problem_operators/time_domain_problem_operator.h
@@ -0,0 +1,28 @@
+#pragma once
+#include "../common/pfem_extras.hpp"
+#include "hephaestus_solvers.h"
+#include "problem_builder_base.h"
+#include "problem_operator_interface.h"
+
+namespace platypus
+{
+
+std::string GetTimeDerivativeName(const std::string & name);
+
+std::vector<std::string> GetTimeDerivativeNames(std::vector<std::string> gridfunction_names);
+
+/// Problem operator for time-dependent problems with no equation system. The user will need to subclass this since the solve is not
+/// implemented.
+class TimeDomainProblemOperator : public mfem::TimeDependentOperator,
+                                  public ProblemOperatorInterface
+{
+public:
+  TimeDomainProblemOperator(platypus::Problem & problem) : ProblemOperatorInterface(problem) {}
+  ~TimeDomainProblemOperator() override = default;
+
+  void SetGridFunctions() override;
+
+  void ImplicitSolve(const double dt, const mfem::Vector & X, mfem::Vector & dX_dt) override {}
+};
+
+} // namespace platypus
\ No newline at end of file
diff --git a/framework/include/mfem/variables/MFEMVariable.h b/framework/include/mfem/variables/MFEMVariable.h
index c997717d..4ae58a1f 100644
--- a/framework/include/mfem/variables/MFEMVariable.h
+++ b/framework/include/mfem/variables/MFEMVariable.h
@@ -1,8 +1,7 @@
 #pragma once
 
 #include "MFEMFESpace.h"
-#include "inputs.hpp"
-#include "gridfunctions.hpp"
+#include "GeneralUserObject.h"
 
 class MFEMVariable : public GeneralUserObject
 {
diff --git a/framework/src/mfem/bcs/MFEMBoundaryCondition.C b/framework/src/mfem/bcs/MFEMBoundaryCondition.C
index 6d596296..75c169bd 100644
--- a/framework/src/mfem/bcs/MFEMBoundaryCondition.C
+++ b/framework/src/mfem/bcs/MFEMBoundaryCondition.C
@@ -27,5 +27,5 @@ MFEMBoundaryCondition::MFEMBoundaryCondition(const InputParameters & parameters)
     bdr_attr[i] = std::stoi(_boundary_names[i]);
   }
   _boundary_condition =
-      std::make_shared<hephaestus::BoundaryCondition>(getParam<std::string>("variable"), bdr_attr);
+      std::make_shared<platypus::BoundaryCondition>(getParam<std::string>("variable"), bdr_attr);
 }
diff --git a/framework/src/mfem/bcs/MFEMScalarDirichletBC.C b/framework/src/mfem/bcs/MFEMScalarDirichletBC.C
index 90d68e61..74aba027 100644
--- a/framework/src/mfem/bcs/MFEMScalarDirichletBC.C
+++ b/framework/src/mfem/bcs/MFEMScalarDirichletBC.C
@@ -15,6 +15,6 @@ MFEMScalarDirichletBC::MFEMScalarDirichletBC(const InputParameters & parameters)
   : MFEMBoundaryCondition(parameters),
     _coef(const_cast<MFEMCoefficient *>(&getUserObject<MFEMCoefficient>("coefficient")))
 {
-  _boundary_condition = std::make_shared<hephaestus::ScalarDirichletBC>(
+  _boundary_condition = std::make_shared<platypus::ScalarDirichletBC>(
       getParam<std::string>("variable"), bdr_attr, _coef->getCoefficient().get());
 }
diff --git a/framework/src/mfem/bcs/MFEMVectorDirichletBC.C b/framework/src/mfem/bcs/MFEMVectorDirichletBC.C
index 496a1aa1..3e806ea7 100644
--- a/framework/src/mfem/bcs/MFEMVectorDirichletBC.C
+++ b/framework/src/mfem/bcs/MFEMVectorDirichletBC.C
@@ -17,6 +17,6 @@ MFEMVectorDirichletBC::MFEMVectorDirichletBC(const InputParameters & parameters)
     _vec_coef(const_cast<MFEMVectorCoefficient *>(
         &getUserObject<MFEMVectorCoefficient>("vector_coefficient")))
 {
-  _boundary_condition = std::make_shared<hephaestus::VectorDirichletBC>(
+  _boundary_condition = std::make_shared<platypus::VectorDirichletBC>(
       getParam<std::string>("variable"), bdr_attr, _vec_coef->getVectorCoefficient().get());
 }
diff --git a/framework/src/mfem/bcs/MFEMVectorNormalIntegratedBC.C b/framework/src/mfem/bcs/MFEMVectorNormalIntegratedBC.C
index 8c4e1fe9..0741e0ef 100644
--- a/framework/src/mfem/bcs/MFEMVectorNormalIntegratedBC.C
+++ b/framework/src/mfem/bcs/MFEMVectorNormalIntegratedBC.C
@@ -19,7 +19,7 @@ MFEMVectorNormalIntegratedBC::MFEMVectorNormalIntegratedBC(const InputParameters
         &getUserObject<MFEMVectorCoefficient>("vector_coefficient")))
 {
 
-  _boundary_condition = std::make_shared<hephaestus::IntegratedBC>(
+  _boundary_condition = std::make_shared<platypus::IntegratedBC>(
       getParam<std::string>("variable"),
       bdr_attr,
       std::make_unique<mfem::BoundaryNormalLFIntegrator>(*_vec_coef->getVectorCoefficient()));
diff --git a/framework/src/mfem/bcs/boundary_conditions.C b/framework/src/mfem/bcs/boundary_conditions.C
new file mode 100644
index 00000000..3260b09b
--- /dev/null
+++ b/framework/src/mfem/bcs/boundary_conditions.C
@@ -0,0 +1,131 @@
+#include "boundary_conditions.h"
+
+namespace platypus
+{
+
+mfem::Array<int>
+BCMap::GetEssentialBdrMarkers(const std::string & name_, mfem::Mesh * mesh_)
+{
+  mfem::Array<int> global_ess_markers(mesh_->bdr_attributes.Max());
+  global_ess_markers = 0;
+  mfem::Array<int> ess_bdrs(mesh_->bdr_attributes.Max());
+  ess_bdrs = 0;
+
+  for (auto const & [name, bc_] : *this)
+  {
+    if (bc_->_name == name_)
+    {
+      auto bc = std::dynamic_pointer_cast<platypus::EssentialBC>(bc_);
+      if (bc != nullptr)
+      {
+        ess_bdrs = bc->GetMarkers(*mesh_);
+        for (auto it = 0; it != mesh_->bdr_attributes.Max(); ++it)
+        {
+          global_ess_markers[it] = std::max(global_ess_markers[it], ess_bdrs[it]);
+        }
+      }
+    }
+  }
+  return global_ess_markers;
+}
+
+void
+BCMap::ApplyEssentialBCs(const std::string & name_,
+                         mfem::Array<int> & ess_tdof_list,
+                         mfem::GridFunction & gridfunc,
+                         mfem::Mesh * mesh_)
+{
+  for (auto const & [name, bc_] : *this)
+  {
+    if (bc_->_name == name_)
+    {
+      auto bc = std::dynamic_pointer_cast<platypus::EssentialBC>(bc_);
+      if (bc != nullptr)
+      {
+        bc->ApplyBC(gridfunc, mesh_);
+      }
+    }
+  }
+  mfem::Array<int> ess_bdr = GetEssentialBdrMarkers(name_, mesh_);
+  gridfunc.FESpace()->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+}
+
+void
+BCMap::ApplyEssentialBCs(const std::string & name_,
+                         mfem::Array<int> & ess_tdof_list,
+                         mfem::ParComplexGridFunction & gridfunc,
+                         mfem::Mesh * mesh_)
+{
+  for (auto const & [name, bc_] : *this)
+  {
+    if (bc_->_name == name_)
+    {
+      auto bc = std::dynamic_pointer_cast<platypus::EssentialBC>(bc_);
+      if (bc != nullptr)
+      {
+        bc->ApplyBC(gridfunc, mesh_);
+      }
+    }
+  }
+  mfem::Array<int> ess_bdr = GetEssentialBdrMarkers(name_, mesh_);
+  gridfunc.FESpace()->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+};
+
+void
+BCMap::ApplyIntegratedBCs(const std::string & name_, mfem::LinearForm & lf, mfem::Mesh * mesh_)
+{
+
+  for (auto const & [name, bc_] : *this)
+  {
+    if (bc_->_name == name_)
+    {
+      auto bc = std::dynamic_pointer_cast<platypus::IntegratedBC>(bc_);
+      if (bc != nullptr)
+      {
+        bc->GetMarkers(*mesh_);
+        bc->ApplyBC(lf);
+      }
+    }
+  }
+};
+
+void
+BCMap::ApplyIntegratedBCs(const std::string & name_,
+                          mfem::ParComplexLinearForm & clf,
+                          mfem::Mesh * mesh_)
+{
+
+  for (auto const & [name, bc_] : *this)
+  {
+    if (bc_->_name == name_)
+    {
+      auto bc = std::dynamic_pointer_cast<platypus::IntegratedBC>(bc_);
+      if (bc != nullptr)
+      {
+        bc->GetMarkers(*mesh_);
+        bc->ApplyBC(clf);
+      }
+    }
+  }
+};
+
+void
+BCMap::ApplyIntegratedBCs(const std::string & name_,
+                          mfem::ParSesquilinearForm & slf,
+                          mfem::Mesh * mesh_)
+{
+  for (auto const & [name, bc_] : *this)
+  {
+    if (bc_->_name == name_)
+    {
+      auto bc = std::dynamic_pointer_cast<platypus::RobinBC>(bc_);
+      if (bc != nullptr)
+      {
+        bc->GetMarkers(*mesh_);
+        bc->ApplyBC(slf);
+      }
+    }
+  }
+};
+
+} // namespace platypus
diff --git a/framework/src/mfem/equation_systems/equation_system.C b/framework/src/mfem/equation_systems/equation_system.C
new file mode 100644
index 00000000..6a3e3366
--- /dev/null
+++ b/framework/src/mfem/equation_systems/equation_system.C
@@ -0,0 +1,431 @@
+#include "equation_system.h"
+
+namespace platypus
+{
+
+EquationSystem::~EquationSystem() { _h_blocks.DeleteAll(); }
+
+bool
+EquationSystem::VectorContainsName(const std::vector<std::string> & the_vector,
+                                   const std::string & name) const
+{
+
+  auto iter = std::find(the_vector.begin(), the_vector.end(), name);
+
+  return (iter != the_vector.end());
+}
+
+void
+EquationSystem::AddTrialVariableNameIfMissing(const std::string & trial_var_name)
+{
+  if (!VectorContainsName(_trial_var_names, trial_var_name))
+  {
+    _trial_var_names.push_back(trial_var_name);
+  }
+}
+
+void
+EquationSystem::AddTestVariableNameIfMissing(const std::string & test_var_name)
+{
+  if (!VectorContainsName(_test_var_names, test_var_name))
+  {
+    _test_var_names.push_back(test_var_name);
+  }
+}
+
+void
+EquationSystem::AddKernel(const std::string & test_var_name,
+                          std::shared_ptr<ParBilinearFormKernel> blf_kernel)
+{
+  AddTestVariableNameIfMissing(test_var_name);
+
+  if (!_blf_kernels_map.Has(test_var_name))
+  {
+    // 1. Create kernels vector.
+    auto kernels = std::make_shared<std::vector<std::shared_ptr<ParBilinearFormKernel>>>();
+
+    // 2. Register with map to prevent leaks.
+    _blf_kernels_map.Register(test_var_name, std::move(kernels));
+  }
+
+  _blf_kernels_map.GetRef(test_var_name).push_back(std::move(blf_kernel));
+}
+
+void
+EquationSystem::AddKernel(const std::string & test_var_name,
+                          std::shared_ptr<ParLinearFormKernel> lf_kernel)
+{
+  AddTestVariableNameIfMissing(test_var_name);
+
+  if (!_lf_kernels_map.Has(test_var_name))
+  {
+    auto kernels = std::make_shared<std::vector<std::shared_ptr<ParLinearFormKernel>>>();
+
+    _lf_kernels_map.Register(test_var_name, std::move(kernels));
+  }
+
+  _lf_kernels_map.GetRef(test_var_name).push_back(std::move(lf_kernel));
+}
+
+void
+EquationSystem::AddKernel(const std::string & test_var_name,
+                          std::shared_ptr<ParNonlinearFormKernel> nlf_kernel)
+{
+  AddTestVariableNameIfMissing(test_var_name);
+
+  if (!_nlf_kernels_map.Has(test_var_name))
+  {
+    auto kernels = std::make_shared<std::vector<std::shared_ptr<ParNonlinearFormKernel>>>();
+
+    _nlf_kernels_map.Register(test_var_name, std::move(kernels));
+  }
+
+  _nlf_kernels_map.GetRef(test_var_name).push_back(std::move(nlf_kernel));
+}
+
+void
+EquationSystem::AddKernel(const std::string & trial_var_name,
+                          const std::string & test_var_name,
+                          std::shared_ptr<ParMixedBilinearFormKernel> mblf_kernel)
+{
+  AddTestVariableNameIfMissing(test_var_name);
+
+  // Register new mblf kernels map if not present for this test variable
+  if (!_mblf_kernels_map_map.Has(test_var_name))
+  {
+    auto kernel_field_map = std::make_shared<
+        platypus::NamedFieldsMap<std::vector<std::shared_ptr<ParMixedBilinearFormKernel>>>>();
+
+    _mblf_kernels_map_map.Register(test_var_name, std::move(kernel_field_map));
+  }
+
+  // Register new mblf kernels map if not present for the test/trial variable
+  // pair
+  if (!_mblf_kernels_map_map.Get(test_var_name)->Has(trial_var_name))
+  {
+    auto kernels = std::make_shared<std::vector<std::shared_ptr<ParMixedBilinearFormKernel>>>();
+
+    _mblf_kernels_map_map.Get(test_var_name)->Register(trial_var_name, std::move(kernels));
+  }
+
+  _mblf_kernels_map_map.GetRef(test_var_name)
+      .Get(trial_var_name)
+      ->push_back(std::move(mblf_kernel));
+}
+
+void
+EquationSystem::ApplyBoundaryConditions(platypus::BCMap & bc_map)
+{
+  _ess_tdof_lists.resize(_test_var_names.size());
+  for (int i = 0; i < _test_var_names.size(); i++)
+  {
+    auto test_var_name = _test_var_names.at(i);
+    // Set default value of gridfunction used in essential BC. Values
+    // overwritten in applyEssentialBCs
+    *(_xs.at(i)) = 0.0;
+    bc_map.ApplyEssentialBCs(
+        test_var_name, _ess_tdof_lists.at(i), *(_xs.at(i)), _test_pfespaces.at(i)->GetParMesh());
+    bc_map.ApplyIntegratedBCs(
+        test_var_name, _lfs.GetRef(test_var_name), _test_pfespaces.at(i)->GetParMesh());
+  }
+}
+void
+EquationSystem::FormLinearSystem(mfem::OperatorHandle & op,
+                                 mfem::BlockVector & trueX,
+                                 mfem::BlockVector & trueRHS)
+{
+
+  // Allocate block operator
+  _h_blocks.DeleteAll();
+  _h_blocks.SetSize(_test_var_names.size(), _test_var_names.size());
+  // Form diagonal blocks.
+  for (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);
+    auto lf = _lfs.Get(test_var_name);
+    mfem::Vector aux_x, aux_rhs;
+    _h_blocks(i, i) = new mfem::HypreParMatrix;
+    blf->FormLinearSystem(
+        _ess_tdof_lists.at(i), *(_xs.at(i)), *lf, *_h_blocks(i, i), aux_x, aux_rhs);
+    trueX.GetBlock(i) = aux_x;
+    trueRHS.GetBlock(i) = aux_rhs;
+  }
+
+  // Form off-diagonal blocks
+  for (int i = 0; i < _test_var_names.size(); i++)
+  {
+    auto test_var_name = _test_var_names.at(i);
+    for (int j = 0; j < _test_var_names.size(); j++)
+    {
+      auto trial_var_name = _test_var_names.at(j);
+
+      mfem::Vector aux_x, aux_rhs;
+      mfem::ParLinearForm aux_lf(_test_pfespaces.at(i));
+      aux_lf = 0.0;
+      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);
+        _h_blocks(i, j) = new mfem::HypreParMatrix;
+        mblf->FormRectangularLinearSystem(_ess_tdof_lists.at(j),
+                                          _ess_tdof_lists.at(i),
+                                          *(_xs.at(j)),
+                                          aux_lf,
+                                          *_h_blocks(i, j),
+                                          aux_x,
+                                          aux_rhs);
+        trueRHS.GetBlock(i) += aux_rhs;
+      }
+    }
+  }
+  // Sync memory
+  for (int i = 0; i < _test_var_names.size(); i++)
+  {
+    trueX.GetBlock(0).SyncAliasMemory(trueX);
+    trueRHS.GetBlock(0).SyncAliasMemory(trueRHS);
+  }
+
+  // Create monolithic matrix
+  op.Reset(mfem::HypreParMatrixFromBlocks(_h_blocks));
+}
+
+void
+EquationSystem::BuildJacobian(mfem::BlockVector & trueX, mfem::BlockVector & trueRHS)
+{
+  height = trueX.Size();
+  width = trueRHS.Size();
+  FormLinearSystem(_jacobian, trueX, trueRHS);
+}
+
+void
+EquationSystem::Mult(const mfem::Vector & x, mfem::Vector & residual) const
+{
+  _jacobian->Mult(x, residual);
+}
+
+mfem::Operator &
+EquationSystem::GetGradient(const mfem::Vector & u) const
+{
+  return *_jacobian;
+}
+
+void
+EquationSystem::RecoverFEMSolution(mfem::BlockVector & trueX,
+                                   platypus::GridFunctions & gridfunctions)
+{
+  for (int i = 0; i < _test_var_names.size(); i++)
+  {
+    auto & test_var_name = _test_var_names.at(i);
+    trueX.GetBlock(i).SyncAliasMemory(trueX);
+    gridfunctions.Get(test_var_name)->Distribute(&(trueX.GetBlock(i)));
+  }
+}
+
+void
+EquationSystem::Init(platypus::GridFunctions & gridfunctions,
+                     const platypus::FESpaces & fespaces,
+                     platypus::BCMap & bc_map,
+                     platypus::Coefficients & coefficients)
+{
+
+  // Add optional kernels to the EquationSystem
+  AddKernels();
+
+  for (auto & test_var_name : _test_var_names)
+  {
+    if (!gridfunctions.Has(test_var_name))
+    {
+      MFEM_ABORT("Test variable " << test_var_name
+                                  << " requested by equation system during initialisation was "
+                                     "not found in gridfunctions");
+    }
+    // Store pointers to variable FESpaces
+    _test_pfespaces.push_back(gridfunctions.Get(test_var_name)->ParFESpace());
+    // Create auxiliary gridfunctions for applying Dirichlet conditions
+    _xs.emplace_back(
+        std::make_unique<mfem::ParGridFunction>(gridfunctions.Get(test_var_name)->ParFESpace()));
+  }
+
+  // Initialise bilinear forms
+
+  for (const auto & [test_var_name, blf_kernels] : _blf_kernels_map)
+  {
+    for (auto & i : *blf_kernels)
+    {
+      i->Init(gridfunctions, fespaces, bc_map, coefficients);
+    }
+  }
+  // Initialise linear form kernels
+  for (const auto & [test_var_name, lf_kernels] : _lf_kernels_map)
+  {
+    for (auto & i : *lf_kernels)
+    {
+      i->Init(gridfunctions, fespaces, bc_map, coefficients);
+    }
+  }
+  // Initialise nonlinear form kernels
+  for (const auto & [test_var_name, nlf_kernels] : _nlf_kernels_map)
+  {
+    for (auto & i : *nlf_kernels)
+    {
+      i->Init(gridfunctions, fespaces, bc_map, coefficients);
+    }
+  }
+  // Initialise mixed bilinear form kernels
+  for (const auto & [test_var_name, mblf_kernels_map] : _mblf_kernels_map_map)
+  {
+    for (const auto & [trial_var_name, mblf_kernels] : *mblf_kernels_map)
+    {
+      for (auto & i : *mblf_kernels)
+      {
+        i->Init(gridfunctions, fespaces, bc_map, coefficients);
+      }
+    }
+  }
+}
+
+void
+EquationSystem::BuildLinearForms(platypus::BCMap & bc_map)
+{
+  // Register linear forms
+  for (int i = 0; i < _test_var_names.size(); i++)
+  {
+    auto test_var_name = _test_var_names.at(i);
+    _lfs.Register(test_var_name, std::make_shared<mfem::ParLinearForm>(_test_pfespaces.at(i)));
+    _lfs.GetRef(test_var_name) = 0.0;
+  }
+  // Apply boundary conditions
+  ApplyBoundaryConditions(bc_map);
+
+  for (auto & test_var_name : _test_var_names)
+  {
+    // Apply kernels
+    auto lf = _lfs.Get(test_var_name);
+    // Assemble. Must be done before applying kernels that add to lf.
+    lf->Assemble();
+
+    if (_lf_kernels_map.Has(test_var_name))
+    {
+      auto lf_kernels = _lf_kernels_map.GetRef(test_var_name);
+
+      for (auto & lf_kernel : lf_kernels)
+      {
+        lf_kernel->Apply(lf);
+      }
+    }
+  }
+}
+
+void
+EquationSystem::BuildBilinearForms()
+{
+  // Register bilinear forms
+  for (int i = 0; i < _test_var_names.size(); i++)
+  {
+    auto test_var_name = _test_var_names.at(i);
+    _blfs.Register(test_var_name, std::make_shared<mfem::ParBilinearForm>(_test_pfespaces.at(i)));
+
+    // Apply kernels
+    auto blf = _blfs.Get(test_var_name);
+    if (_blf_kernels_map.Has(test_var_name))
+    {
+      auto blf_kernels = _blf_kernels_map.GetRef(test_var_name);
+
+      for (auto & blf_kernel : blf_kernels)
+      {
+        blf_kernel->Apply(blf);
+      }
+    }
+    // Assemble
+    blf->Assemble();
+  }
+}
+
+void
+EquationSystem::BuildMixedBilinearForms()
+{
+  // Register mixed linear forms. Note that not all combinations may
+  // have a kernel
+
+  // Create mblf for each test/trial pair
+  for (int i = 0; i < _test_var_names.size(); i++)
+  {
+    auto test_var_name = _test_var_names.at(i);
+    auto test_mblfs = std::make_shared<platypus::NamedFieldsMap<mfem::ParMixedBilinearForm>>();
+    for (int j = 0; j < _test_var_names.size(); j++)
+    {
+      auto trial_var_name = _test_var_names.at(j);
+
+      // Register MixedBilinearForm if kernels exist for it, and assemble
+      // kernels
+      if (_mblf_kernels_map_map.Has(test_var_name) &&
+          _mblf_kernels_map_map.Get(test_var_name)->Has(trial_var_name))
+      {
+        auto mblf_kernels = _mblf_kernels_map_map.GetRef(test_var_name).GetRef(trial_var_name);
+        auto mblf = std::make_shared<mfem::ParMixedBilinearForm>(_test_pfespaces.at(j),
+                                                                 _test_pfespaces.at(i));
+        // Apply all mixed kernels with this test/trial pair
+        for (auto & mblf_kernel : mblf_kernels)
+        {
+          mblf_kernel->Apply(mblf.get());
+        }
+        // Assemble mixed bilinear forms
+        mblf->Assemble();
+        // Register mixed bilinear forms associated with a single trial variable
+        // for the current test variable
+        test_mblfs->Register(trial_var_name, mblf);
+      }
+    }
+    // Register all mixed bilinear form sets associated with a single test
+    // variable
+    _mblfs.Register(test_var_name, test_mblfs);
+  }
+}
+
+void
+EquationSystem::BuildEquationSystem(platypus::BCMap & bc_map)
+{
+  BuildLinearForms(bc_map);
+  BuildBilinearForms();
+  BuildMixedBilinearForms();
+}
+
+TimeDependentEquationSystem::TimeDependentEquationSystem() : _dt_coef(1.0) {}
+
+void
+TimeDependentEquationSystem::AddTrialVariableNameIfMissing(const std::string & var_name)
+{
+  EquationSystem::AddTrialVariableNameIfMissing(var_name);
+  std::string var_time_derivative_name = GetTimeDerivativeName(var_name);
+  if (std::find(_trial_var_time_derivative_names.begin(),
+                _trial_var_time_derivative_names.end(),
+                var_time_derivative_name) == _trial_var_time_derivative_names.end())
+  {
+    _trial_var_time_derivative_names.push_back(var_time_derivative_name);
+  }
+}
+
+void
+TimeDependentEquationSystem::SetTimeStep(double dt)
+{
+  if (fabs(dt - _dt_coef.constant) > 1.0e-12 * dt)
+  {
+    _dt_coef.constant = dt;
+    for (auto test_var_name : _test_var_names)
+    {
+      auto blf = _blfs.Get(test_var_name);
+      blf->Update();
+      blf->Assemble();
+    }
+  }
+}
+
+void
+TimeDependentEquationSystem::UpdateEquationSystem(platypus::BCMap & bc_map)
+{
+  BuildLinearForms(bc_map);
+  BuildBilinearForms();
+  BuildMixedBilinearForms();
+}
+
+} // namespace platypus
diff --git a/framework/src/mfem/kernels/MFEMDiffusionKernel.C b/framework/src/mfem/kernels/MFEMDiffusionKernel.C
index bdd2ee3f..0c22086d 100644
--- a/framework/src/mfem/kernels/MFEMDiffusionKernel.C
+++ b/framework/src/mfem/kernels/MFEMDiffusionKernel.C
@@ -20,6 +20,6 @@ MFEMDiffusionKernel::MFEMDiffusionKernel(const InputParameters & parameters)
   : MFEMBilinearFormKernel(parameters),
     _kernel_params{{{"VariableName", getParam<std::string>("variable")},
                     {"CoefficientName", getParam<std::string>("coefficient")}}},
-    _kernel{std::make_shared<hephaestus::DiffusionKernel>(_kernel_params)}
+    _kernel{std::make_shared<platypus::DiffusionKernel>(_kernel_params)}
 {
 }
diff --git a/framework/src/mfem/problem/MFEMProblem.C b/framework/src/mfem/problem/MFEMProblem.C
index 5de8121f..66dd591a 100644
--- a/framework/src/mfem/problem/MFEMProblem.C
+++ b/framework/src/mfem/problem/MFEMProblem.C
@@ -23,7 +23,6 @@ MFEMProblem::MFEMProblem(const InputParameters & params)
     _outputs(),
     _exec_params()
 {
-  hephaestus::logger.set_level(spdlog::level::info);
 }
 
 MFEMProblem::~MFEMProblem() {}
@@ -69,13 +68,13 @@ MFEMProblem::initialSetup()
   // NB: set to false to avoid reconstructing problem operator.
   mfem_problem_builder->FinalizeProblem(false);
 
-  hephaestus::InputParameters exec_params;
+  platypus::InputParameters exec_params;
 
   Transient * _moose_executioner = dynamic_cast<Transient *>(_app.getExecutioner());
   if (_moose_executioner != nullptr)
   {
     auto mfem_transient_problem_builder =
-        std::dynamic_pointer_cast<hephaestus::TimeDomainProblemBuilder>(mfem_problem_builder);
+        std::dynamic_pointer_cast<platypus::TimeDomainProblemBuilder>(mfem_problem_builder);
     if (mfem_transient_problem_builder == nullptr)
     {
       mooseError("Specified formulation does not support Transient executioners");
@@ -87,15 +86,14 @@ MFEMProblem::initialSetup()
     exec_params.SetParam("TimeStep", float(dt()));
     exec_params.SetParam("EndTime", float(_moose_executioner->endTime()));
     exec_params.SetParam("VisualisationSteps", getParam<int>("vis_steps"));
-    exec_params.SetParam("Problem",
-                         static_cast<hephaestus::TimeDomainProblem *>(mfem_problem.get()));
+    exec_params.SetParam("Problem", static_cast<platypus::TimeDomainProblem *>(mfem_problem.get()));
 
-    executioner = std::make_unique<hephaestus::TransientExecutioner>(exec_params);
+    executioner = std::make_unique<platypus::TransientExecutioner>(exec_params);
   }
   else if (dynamic_cast<Steady *>(_app.getExecutioner()))
   {
     auto mfem_steady_problem_builder =
-        std::dynamic_pointer_cast<hephaestus::SteadyStateProblemBuilder>(mfem_problem_builder);
+        std::dynamic_pointer_cast<platypus::SteadyStateProblemBuilder>(mfem_problem_builder);
     if (mfem_steady_problem_builder == nullptr)
     {
       mooseError("Specified formulation does not support Steady executioners");
@@ -104,9 +102,9 @@ MFEMProblem::initialSetup()
     mfem_problem = mfem_steady_problem_builder->ReturnProblem();
 
     exec_params.SetParam("Problem",
-                         static_cast<hephaestus::SteadyStateProblem *>(mfem_problem.get()));
+                         static_cast<platypus::SteadyStateProblem *>(mfem_problem.get()));
 
-    executioner = std::make_unique<hephaestus::SteadyExecutioner>(exec_params);
+    executioner = std::make_unique<platypus::SteadyExecutioner>(exec_params);
   }
   else
   {
@@ -130,7 +128,7 @@ MFEMProblem::externalSolve()
     return;
   }
 
-  auto * transient_mfem_exec = dynamic_cast<hephaestus::TransientExecutioner *>(executioner.get());
+  auto * transient_mfem_exec = dynamic_cast<platypus::TransientExecutioner *>(executioner.get());
   if (transient_mfem_exec != nullptr)
   {
     transient_mfem_exec->_t_step = dt();
@@ -174,7 +172,7 @@ MFEMProblem::addMaterial(const std::string & kernel_name,
   for (unsigned int bid = 0; bid < mfem_material.blocks.size(); ++bid)
   {
     int block = std::stoi(mfem_material.blocks[bid]);
-    hephaestus::Subdomain mfem_subdomain(name, block);
+    platypus::Subdomain mfem_subdomain(name, block);
     mfem_material.storeCoefficients(mfem_subdomain);
     _coefficients._subdomains.push_back(mfem_subdomain);
   }
@@ -188,14 +186,6 @@ MFEMProblem::addCoefficient(const std::string & user_object_name,
   FEProblemBase::addUserObject(user_object_name, name, parameters);
   MFEMCoefficient * mfem_coef(&getUserObject<MFEMCoefficient>(name));
   _coefficients._scalars.Register(name, mfem_coef->getCoefficient());
-
-  // Add associated auxsolvers for CoupledCoefficients
-  auto coupled_coef = std::dynamic_pointer_cast<hephaestus::CoupledCoefficient>(
-      _coefficients._scalars.GetShared(name));
-  if (coupled_coef != nullptr)
-  {
-    mfem_problem_builder->AddAuxSolver(name, std::move(coupled_coef));
-  }
 }
 
 void
@@ -273,16 +263,8 @@ MFEMProblem::addAuxKernel(const std::string & kernel_name,
 {
   std::string base_auxkernel = parameters.get<std::string>("_moose_base");
 
-  if (base_auxkernel == "MFEMAuxKernel") // MFEM auxsolver.
-  {
-    FEProblemBase::addUserObject(kernel_name, name, parameters);
-    MFEMAuxSolver * mfem_auxsolver(&getUserObject<MFEMAuxSolver>(name));
-
-    mfem_problem_builder->AddPostprocessor(name, mfem_auxsolver->getAuxSolver());
-    mfem_auxsolver->storeCoefficients(_coefficients);
-  }
-  else if (base_auxkernel == "AuxKernel" || base_auxkernel == "VectorAuxKernel" ||
-           base_auxkernel == "ArrayAuxKernel") // MOOSE auxkernels.
+  if (base_auxkernel == "AuxKernel" || base_auxkernel == "VectorAuxKernel" ||
+      base_auxkernel == "ArrayAuxKernel") // MOOSE auxkernels.
   {
     FEProblemBase::addAuxKernel(kernel_name, name, parameters);
   }
diff --git a/framework/src/mfem/problem_operators/equation_system_problem_operator.C b/framework/src/mfem/problem_operators/equation_system_problem_operator.C
new file mode 100644
index 00000000..bc9c48e0
--- /dev/null
+++ b/framework/src/mfem/problem_operators/equation_system_problem_operator.C
@@ -0,0 +1,20 @@
+#include "equation_system_problem_operator.h"
+
+namespace platypus
+{
+void
+EquationSystemProblemOperator::SetGridFunctions()
+{
+  _trial_var_names = GetEquationSystem()->_trial_var_names;
+  ProblemOperator::SetGridFunctions();
+}
+
+void
+EquationSystemProblemOperator::Init(mfem::Vector & X)
+{
+  ProblemOperator::Init(X);
+
+  GetEquationSystem()->BuildEquationSystem(_problem._bc_map);
+}
+
+} // namespace platypus
\ No newline at end of file
diff --git a/framework/src/mfem/problem_operators/problem_operator.C b/framework/src/mfem/problem_operators/problem_operator.C
new file mode 100644
index 00000000..e6f3c034
--- /dev/null
+++ b/framework/src/mfem/problem_operators/problem_operator.C
@@ -0,0 +1,13 @@
+#include "problem_operator.h"
+
+namespace platypus
+{
+
+void
+ProblemOperator::SetGridFunctions()
+{
+  ProblemOperatorInterface::SetGridFunctions();
+  width = height = _true_offsets[_trial_variables.size()];
+};
+
+} // namespace platypus
\ No newline at end of file
diff --git a/framework/src/mfem/problem_operators/problem_operator_interface.C b/framework/src/mfem/problem_operators/problem_operator_interface.C
new file mode 100644
index 00000000..170b4f2c
--- /dev/null
+++ b/framework/src/mfem/problem_operators/problem_operator_interface.C
@@ -0,0 +1,43 @@
+#include "problem_operator_interface.h"
+
+namespace platypus
+{
+void
+ProblemOperatorInterface::SetGridFunctions()
+{
+  _trial_variables = _problem._gridfunctions.Get(_trial_var_names);
+
+  // Set operator size and block structure
+  _block_true_offsets.SetSize(_trial_variables.size() + 1);
+  _block_true_offsets[0] = 0;
+  for (unsigned int ind = 0; ind < _trial_variables.size(); ++ind)
+  {
+    _block_true_offsets[ind + 1] = _trial_variables.at(ind)->ParFESpace()->TrueVSize();
+  }
+  _block_true_offsets.PartialSum();
+
+  _true_offsets.SetSize(_trial_variables.size() + 1);
+  _true_offsets[0] = 0;
+  for (unsigned int ind = 0; ind < _trial_variables.size(); ++ind)
+  {
+    _true_offsets[ind + 1] = _trial_variables.at(ind)->ParFESpace()->GetVSize();
+  }
+  _true_offsets.PartialSum();
+
+  _true_x.Update(_block_true_offsets);
+  _true_rhs.Update(_block_true_offsets);
+}
+
+void
+ProblemOperatorInterface::Init(mfem::Vector & X)
+{
+  for (size_t i = 0; i < _trial_variables.size(); ++i)
+  {
+    mfem::ParGridFunction * trial_var = _trial_variables.at(i);
+
+    trial_var->MakeRef(trial_var->ParFESpace(), const_cast<mfem::Vector &>(X), _true_offsets[i]);
+    *trial_var = 0.0;
+  }
+}
+
+}
\ No newline at end of file
diff --git a/framework/src/mfem/problem_operators/time_domain_equation_system_problem_operator.C b/framework/src/mfem/problem_operators/time_domain_equation_system_problem_operator.C
new file mode 100644
index 00000000..099fb086
--- /dev/null
+++ b/framework/src/mfem/problem_operators/time_domain_equation_system_problem_operator.C
@@ -0,0 +1,61 @@
+#include "time_domain_equation_system_problem_operator.h"
+
+namespace platypus
+{
+
+void
+TimeDomainEquationSystemProblemOperator::SetGridFunctions()
+{
+  _trial_var_names = GetEquationSystem()->_trial_var_names;
+  _trial_variable_time_derivatives =
+      _problem._gridfunctions.Get(GetEquationSystem()->_trial_var_time_derivative_names);
+
+  TimeDomainProblemOperator::SetGridFunctions();
+}
+
+void
+TimeDomainEquationSystemProblemOperator::Init(mfem::Vector & X)
+{
+  TimeDomainProblemOperator::Init(X);
+
+  // Define material property coefficients
+  for (size_t i = 0; i < _trial_variables.size(); ++i)
+  {
+    *(_trial_variable_time_derivatives.at(i)) = 0.0;
+  }
+
+  GetEquationSystem()->BuildEquationSystem(_problem._bc_map);
+}
+
+void
+TimeDomainEquationSystemProblemOperator::ImplicitSolve(const double dt,
+                                                       const mfem::Vector & X,
+                                                       mfem::Vector & dX_dt)
+{
+  dX_dt = 0.0;
+  for (unsigned int ind = 0; ind < _trial_variables.size(); ++ind)
+  {
+    _trial_variables.at(ind)->MakeRef(
+        _trial_variables.at(ind)->ParFESpace(), const_cast<mfem::Vector &>(X), _true_offsets[ind]);
+    _trial_variable_time_derivatives.at(ind)->MakeRef(
+        _trial_variable_time_derivatives.at(ind)->ParFESpace(), dX_dt, _true_offsets[ind]);
+  }
+  _problem._coefficients.SetTime(GetTime());
+  BuildEquationSystemOperator(dt);
+
+  _problem._nonlinear_solver->SetSolver(*_problem._jacobian_solver);
+  _problem._nonlinear_solver->SetOperator(*GetEquationSystem());
+  _problem._nonlinear_solver->Mult(_true_rhs, _true_x);
+
+  GetEquationSystem()->RecoverFEMSolution(_true_x, _problem._gridfunctions);
+}
+
+void
+TimeDomainEquationSystemProblemOperator::BuildEquationSystemOperator(double dt)
+{
+  GetEquationSystem()->SetTimeStep(dt);
+  GetEquationSystem()->UpdateEquationSystem(_problem._bc_map);
+  GetEquationSystem()->BuildJacobian(_true_x, _true_rhs);
+}
+
+} // namespace platypus
\ No newline at end of file
diff --git a/framework/src/mfem/problem_operators/time_domain_problem_operator.C b/framework/src/mfem/problem_operators/time_domain_problem_operator.C
new file mode 100644
index 00000000..68cbf7ae
--- /dev/null
+++ b/framework/src/mfem/problem_operators/time_domain_problem_operator.C
@@ -0,0 +1,30 @@
+#include "time_domain_problem_operator.h"
+
+namespace platypus
+{
+
+std::string
+GetTimeDerivativeName(const std::string & name)
+{
+  return std::string("d") + name + std::string("_dt");
+}
+
+std::vector<std::string>
+GetTimeDerivativeNames(std::vector<std::string> gridfunction_names)
+{
+  std::vector<std::string> time_derivative_names;
+  for (auto & gridfunction_name : gridfunction_names)
+  {
+    time_derivative_names.push_back(GetTimeDerivativeName(gridfunction_name));
+  }
+  return time_derivative_names;
+}
+
+void
+TimeDomainProblemOperator::SetGridFunctions()
+{
+  ProblemOperatorInterface::SetGridFunctions();
+  width = height = _true_offsets[_trial_variables.size()];
+}
+
+} // namespace platypus
\ No newline at end of file