diff --git a/Field.hpp b/Field.hpp
index 5e97e51bf20f02fb6c858133c0879cf4ed29f8ba..a5e900759aafa825ad2f2b7e85a188d03c4cab06 100644
--- a/Field.hpp
+++ b/Field.hpp
@@ -99,6 +99,39 @@ public:
         return *this;
     }
 
+  /**
+   * @brief Compare this field to another field, with absolute
+   * tolerance tol.
+   * @return bool
+   */
+  bool compare(Field<TType, TState> &rhs, double tol)
+  {
+    const std::vector<BlockAttributes>& rhs_blocks = rhs.GetBlocks();
+    TType *store = GetStorage().GetCPUPtr();
+    TType *rhs_store = rhs.GetStorage().GetCPUPtr();
+    if (rhs_blocks.size() != block_attributes.size()) return false;
+
+    size_t i {0};
+    for (size_t bl = 0; bl < block_attributes.size(); ++bl){
+
+      size_t num_elements = block_attributes[bl].num_elements;
+      size_t num_pts = block_attributes[bl].num_pts;
+
+      // Check that each block have the same structure
+      if (num_elements != rhs_blocks[bl].num_elements) return false;
+      if (num_pts != rhs_blocks[bl].num_pts) return false;
+
+      // Compare elements in the blocks
+      for (size_t el = 0 ; el < num_elements; ++el) {
+	for (size_t coeff = 0; coeff < num_pts; ++coeff) {
+	  i = coeff + el * num_pts + bl * num_pts * num_elements;
+	  if (std::abs(store[i] - rhs_store[i]) > tol) return false;
+	}
+      }
+    }
+    return true;
+  }
+
     /**
      * @brief Static templated creation method.
      *
diff --git a/tests/init_fields.hpp b/tests/init_fields.hpp
index 9aeb7f82f39bc72bb5219feff10ef56e3139e456..d623d935d7584e4f50323e4cc454c497358d0748 100644
--- a/tests/init_fields.hpp
+++ b/tests/init_fields.hpp
@@ -32,6 +32,7 @@ class InitFields
 public:
     Field<double, stateIn> *fixt_in;
     Field<double, stateOut> *fixt_out;
+    Field<double, stateOut> *fixt_expected;
     MultiRegions::ExpListSharedPtr fixt_explist{nullptr};
     ~InitFields()
     {
@@ -67,8 +68,10 @@ public:
         // Create two Field objects with a MemoryRegionCPU backend by default
         auto f_in  = Field<double, stateIn>::create(blocks_in);
         auto f_out = Field<double, stateOut>::create(blocks_out);
+	auto f_expected = Field<double, stateOut>::create(blocks_out);
         fixt_in    = new Field<double, stateIn>(std::move(f_in));
         fixt_out   = new Field<double, stateOut>(std::move(f_out));
+	fixt_expected = new Field<double, stateOut>(std::move(f_expected));
     }
 
 protected:
diff --git a/tests/test_bwdtrans.cpp b/tests/test_bwdtrans.cpp
index 5116addecf30f5e78ce15aa7b77a1f4725423c6e..70c2a51f03456facac226ff46a19b869e4820c8e 100644
--- a/tests/test_bwdtrans.cpp
+++ b/tests/test_bwdtrans.cpp
@@ -63,8 +63,28 @@ BOOST_FIXTURE_TEST_CASE(bwdtrans, Line)
         }
     }
 
-    BwdTrans<>::create(fixt_explist, "StdMat")->apply(*fixt_in, *fixt_out);
+    // TODO: Initialise expected solution
+    x = fixt_expected->GetStorage().GetCPUPtr();
 
-    double *y = fixt_out->GetStorage().GetCPUPtr();
-    BOOST_TEST(y[0] == 1.0);
+    // For each element, initialise first coefficient to zero and rest
+    // to 1.
+    for (auto const &block : fixt_expected->GetBlocks())
+    {
+        for (size_t el = 0; el < block.num_elements; ++el)
+        {
+            for (size_t coeff = 0; coeff < block.num_pts; ++coeff)
+            {
+                if (coeff == 0) {
+                    *(x++) = 1.0;
+                }
+                else {
+                    *(x++) = 0.0;
+                }
+            }
+        }
+    }
+
+    BwdTrans<>::create(fixt_explist, "StdMat")->apply(*fixt_in, *fixt_out);
+    double TOL  {0.01};
+    BOOST_TEST( fixt_out->compare(*fixt_expected, TOL));
 }