diff --git a/include/problem/MFEMProblem.h b/include/problem/MFEMProblem.h
index 07156557..5a60f028 100644
--- a/include/problem/MFEMProblem.h
+++ b/include/problem/MFEMProblem.h
@@ -126,6 +126,10 @@ class MFEMProblem : public ExternalProblem
                    const std::string & name,
                    InputParameters & parameters) override;
 
+  void addInitialCondition(const std::string & ic_name,
+                           const std::string & name,
+                           InputParameters & parameters) override;
+
   /**
    * Method called in AddMFEMPreconditionerAction which will create the solver.
    */
@@ -219,6 +223,16 @@ class MFEMProblem : public ExternalProblem
   std::optional<std::reference_wrapper<mfem::ParGridFunction const>>
   getMeshDisplacementGridFunction();
 
+  /**
+   * @returns a shared pointer to an MFEM coefficient created using the [Coefficients] syntax
+   */
+  std::shared_ptr<mfem::Coefficient> getCoefficient(const std::string & name);
+
+  /**
+   * @returns a shared pointer to an MFEM parallel grid function
+   */
+  std::shared_ptr<mfem::ParGridFunction> getGridFunction(const std::string & name);
+
 protected:
   /**
    * Template method for adding kernels. We can only add kernels using equation system problem
diff --git a/src/base/PlatypusApp.C b/src/base/PlatypusApp.C
index 0bd90620..ac7da8ff 100644
--- a/src/base/PlatypusApp.C
+++ b/src/base/PlatypusApp.C
@@ -38,7 +38,8 @@ associateSyntaxInner(Syntax & syntax, ActionFactory & /*action_factory*/)
   addTaskDependency("add_material", "add_mfem_coefficients");
   addTaskDependency("add_mfem_coefficients", "add_variable");
   addTaskDependency("add_mfem_coefficients", "add_aux_variable");
-  addTaskDependency("add_mfem_coefficients", "add_ic");
+  addTaskDependency("add_ic", "add_mfem_coefficients");
+  addTaskDependency("add_mfem_coefficients", "add_function");
 
   // add vector coefficients
   registerMooseObjectTask("add_mfem_vector_coefficients", MFEMVectorCoefficient, false);
diff --git a/src/problem/MFEMProblem.C b/src/problem/MFEMProblem.C
index b225df27..98a5c3ee 100644
--- a/src/problem/MFEMProblem.C
+++ b/src/problem/MFEMProblem.C
@@ -1,4 +1,5 @@
 #include "MFEMProblem.h"
+#include "MFEMScalarIC.h"
 
 #include <vector>
 #include <algorithm>
@@ -421,3 +422,24 @@ MFEMProblem::mesh()
               "Please choose the MFEMMesh mesh type for an MFEMProblem\n");
   return (MFEMMesh &)_mesh;
 }
+
+std::shared_ptr<mfem::Coefficient>
+MFEMProblem::getCoefficient(const std::string & name)
+{
+  return getUserObject<MFEMCoefficient>(name).getCoefficient();
+}
+
+std::shared_ptr<mfem::ParGridFunction>
+MFEMProblem::getGridFunction(const std::string & name)
+{
+  return getUserObject<MFEMVariable>(name).getGridFunction();
+}
+
+void
+MFEMProblem::addInitialCondition(const std::string & ic_name,
+                                 const std::string & name,
+                                 InputParameters & parameters)
+{
+  FEProblemBase::addUserObject(ic_name, name, parameters);
+  getUserObject<MFEMScalarIC>(name); // error check
+}
diff --git a/test/tests/kernels/diffusion.i b/test/tests/kernels/diffusion.i
index 188de73f..c1aa4034 100644
--- a/test/tests/kernels/diffusion.i
+++ b/test/tests/kernels/diffusion.i
@@ -23,12 +23,20 @@
   []
 []
 
+[ICs]
+  [diffused_ic]
+    type = MFEMScalarIC
+    coefficient = one
+    variable = diffused
+  []
+[]
+
 [Functions]
-  [value_bottom]
+  [one]
     type = ParsedFunction
     expression = 1.0
   []
-  [value_top]
+  [zero]
     type = ParsedFunction
     expression = 0.0
   []
@@ -39,13 +47,13 @@
     type = MFEMScalarDirichletBC
     variable = diffused
     boundary = '1'
-    coefficient = BottomValue
+    coefficient = one
   []
   [low_terminal]
     type = MFEMScalarDirichletBC
     variable = diffused
     boundary = '2'
-    coefficient = TopValue
+    coefficient = zero
   []
 []
 
@@ -58,13 +66,13 @@
 []
 
 [Coefficients]
-  [TopValue]
+  [zero]
     type = MFEMFunctionCoefficient
-    function = value_top
+    function = zero
   []
-  [BottomValue]
+  [one]
     type = MFEMFunctionCoefficient
-    function = value_bottom
+    function = one
   []
 []
 
@@ -86,7 +94,7 @@
   type = MFEMHypreGMRES
   preconditioner = boomeramg
   l_tol = 1e-16
-  l_max_its = 1000  
+  l_max_its = 1000
 []
 
 [Executioner]