diff --git a/examples/diffusion/mesh1D.py b/examples/diffusion/mesh1D.py
index cc5adac08a..6c64b8a65d 100755
--- a/examples/diffusion/mesh1D.py
+++ b/examples/diffusion/mesh1D.py
@@ -731,7 +731,7 @@
 ...     input("No-flux - steady-state failure. \
 ... Press <return> to proceed...")
 
->>> print(numerix.allclose(phi, 0.2, atol=1e-5))
+>>> print(numerix.allclose(phi, 0.2, atol=1e-5)) # doctest: +NOT_TRILINOS_SOLVER
 False
 
 .. image:: /figures/examples/diffusion/mesh1D-noflux_steady_fail.*
@@ -742,6 +742,8 @@
 Depending on the solver, we find that the value may be uniformly zero,
 infinity, or NaN, or the solver may just fail!
 What happened to our no-flux boundary conditions?
+Trilinos actually manages to get the correct solution, but this should not
+be relied on; this problem has an infinite number of solutions.
 
 The problem is that in the implicit discretization of :math:`\nabla\cdot(D\nabla\phi) = 0`,
 
diff --git a/fipy/solvers/__init__.py b/fipy/solvers/__init__.py
index 01180b9267..613e68b9ec 100644
--- a/fipy/solvers/__init__.py
+++ b/fipy/solvers/__init__.py
@@ -190,6 +190,12 @@ def _import_mesh_matrices(suite):
                  test=lambda: solver_suite != 'pyamgx',
                  why="the PyAMGX solver is being used.",
                  skipWarning=True)
+
+register_skipper(flag='NOT_TRILINOS_SOLVER',
+                 test=lambda: solver_suite not in ['trilinos', 'no-pysparse'],
+                 why="the Trilinos solvers are being used.",
+                 skipWarning=True)
+
 del register_skipper
 
 _log.info("Solver suite is %s", solver_suite)