diff --git a/.gitignore b/.gitignore
index eb4b515..f9ec7e0 100644
--- a/.gitignore
+++ b/.gitignore
@@ -40,5 +40,9 @@ manual/
 # Emacs configuration
 .dir-locals.el
 
+# Output
+*.dat
+*.png
+
 # Stupid MacOS
-.DS_Store
\ No newline at end of file
+.DS_Store
diff --git a/CHANGELOG.md b/CHANGELOG.md
index c180acf..fbe858c 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -6,6 +6,7 @@ This is the change log for `NuCFD`, it is based on the format suggested by [Keep
 
 ### Added
 
+- Added example code to compute derivative of a function [010989b].
 - Added `FORD` documentation system [c780fe1].
 - Added tridiagonal solver with support for periodic problems [f5d254f].
 - Added simple library to support writing tests [a5814c8].
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 2f62c41..d0d8dc8 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -50,8 +50,14 @@ endif()
 add_subdirectory(src)
 
 ## Testing
-enable_testing()
-add_subdirectory(tests)
+# option(BUILD_TESTING "Build the testing tree" OFF) # Disables testing by default
+include(CTest)
+if(BUILD_TESTING)
+  add_subdirectory(tests)
+endif()
+
+## Examples (optional)
+add_subdirectory(examples EXCLUDE_FROM_ALL)
 
 ## Documentation
 add_custom_target(doc ford
diff --git a/NuCFD.md b/NuCFD.md
index 0171cd1..dfe4324 100644
--- a/NuCFD.md
+++ b/NuCFD.md
@@ -4,6 +4,7 @@ author: Paul Bartholomew
 source: true
 src_dir: ./src
          ./tests
+		 ./examples
 output_dir: ./manual
 page_dir: ./docs
 ---
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
new file mode 100644
index 0000000..330c1f8
--- /dev/null
+++ b/examples/CMakeLists.txt
@@ -0,0 +1,18 @@
+## examples/CMakeLists.txt
+#
+## Description
+#
+# CMakeLists file for building the NuCFD project examples.
+#
+## LICENSE
+#
+# SPDX-License-Identifier: BSD-3-Clause
+#
+
+project(examples)
+
+add_executable(compute_derivatives compute-derivatives.f90)
+target_link_libraries(compute_derivatives nucfd)
+
+add_custom_target(${PROJECT_NAME})
+add_dependencies(${PROJECT_NAME} compute_derivatives)
diff --git a/examples/compute-derivatives.f90 b/examples/compute-derivatives.f90
new file mode 100644
index 0000000..f22b77d
--- /dev/null
+++ b/examples/compute-derivatives.f90
@@ -0,0 +1,219 @@
+! examples/compute-derivatives.f90
+!
+!! Example program to compute derivatives.
+!
+! SPDX-License-Identifier: BSD-3-Clause
+
+program compute_derivatives
+  !! A program that uses non-uniform compact finite difference schemes to compute derivatives on
+  !! (non-)uniform meshes.
+
+  use nucfd_types
+  use nucfd_deriv
+  use nucfd_trid_solver
+  
+  implicit none
+
+  real, parameter :: pi = 4.0 * atan(1.0)
+  real, parameter :: alpha = 1.0 / 3.0
+  
+  integer :: n
+  real :: L
+  integer, dimension(:), allocatable :: block_spacings
+  logical :: periodic
+  real, dimension(:), allocatable :: x ! 1D "grid"
+
+  type(nucfd_index_stencil) :: stencil
+  integer :: width, centre
+
+  integer :: idx, i
+
+  real, dimension(:), allocatable :: f
+  real, dimension(:), allocatable :: dfdx
+
+  real, dimension(:), allocatable :: a, b, c, r
+  
+  n = 66
+  L = 1.0
+  block_spacings = (/ 1, 2, 1 /)
+  periodic = .true.
+  call build_grid(n, L, block_spacings, periodic, x)
+
+  print *, "+++ Built grid +++"
+  print *, "- lbounds: ", lbound(x)
+  print *, "- ubounds: ", ubound(x)
+  print *, "- 1st/last node: ", x(1), x(n)
+  if (periodic) then
+     print *, "- periodic node: ", x(n + 1)
+  end if
+
+  allocate(f, mold=x)
+  f(:) = sin((x(:) / L) * (2.0 * pi))
+  
+  allocate(a(n))
+  allocate(b(n))
+  allocate(c(n))
+  allocate(r(n))
+  allocate(dfdx(n))
+  
+  width = 5
+  centre = 3
+  call create_stencil(width, centre, stencil)
+
+  print *, "+++ Computing RHS +++"
+  select type (indices => stencil%stencil)
+  type is(integer)
+     do idx = 1, n
+        ! Move stencil
+        do i = lbound(indices, 1), ubound(indices, 1)
+           indices(i) = idx + i 
+        end do
+        
+        call deriv_rhs(f, stencil, x, r(idx))
+     end do
+  end select
+
+  print *, "+++ Solving system +++"
+  a(:) = alpha
+  b(:) = 1.0
+  c(:) = alpha
+  call solve_cyclic(a, b, c, r, dfdx)
+
+  open (unit=10, file="dfdx.dat")
+  do i = 1, n
+     write(10, *) x(i), (2.0 * pi / L) * cos((x(i) / L) * (2.0 * pi)), dfdx(i)
+  end do
+  close(10)
+
+  deallocate(f, dfdx)
+  deallocate(x)
+
+  deallocate(a, b, c, r)
+  
+contains
+
+  subroutine build_grid(n, L, spacings, periodic, x)
+    !! Crude subroutine to build a 1-D grid.
+    !! The grid is described by the number of nodes, length and an array of relative grid spacings
+    !! for each block (constant within a block), returning a 1-D array of node centres.
+    !
+    !! @note The number of cells (n - 1) must be divisible by the number of blocks specified by the
+    !! spacings array otherwise an error will be raised.
+    
+    integer, intent(in) :: n                          !! The number of nodes in the mesh
+    real, intent(in) :: L                             !! The length of the domain
+    integer, dimension(:), intent(in) :: spacings     !! An array of the relative (integer) spacings
+                                                      !! for each block
+    logical, intent(in) :: periodic                   !! Flag indicating periodicity of mesh
+    real, dimension(:), allocatable, intent(out) :: x !! The grid node locations
+
+    integer :: ncells
+    integer :: nblocks
+    integer :: block_size
+    real :: block_length
+    real :: h
+
+    integer :: i
+    integer :: b
+    integer :: idx
+
+    integer, parameter :: lghost = 2 ! Ghost node count at left boundary
+    integer, parameter :: rghost = 2 ! Ghost node count at right boundary
+
+    allocate(x(1-lghost:n+rghost))
+
+    if (.not. periodic) then
+       ncells = n - 1
+    else
+       ncells = n
+    end if
+    nblocks = size(spacings, 1)
+    if (.not. periodic) then
+       block_size = ncells / nblocks
+    else
+       block_size = ncells / nblocks
+    end if
+    if (.not. periodic) then
+       if (block_size * nblocks + 1 /= n) then
+          print *, "Error: invalid grid specification!"
+          print *, " - Grid size n should be specified so that n - 1 is divisible by number of blocks"
+          print *, "   e.g. for 2 blocks a grid size of 11 is valid, 10 is not."
+          error stop
+       end if
+    else
+       if (block_size * nblocks /= n) then
+          print *, "Error: invalid grid specification!"
+          error stop
+       end if
+    end if
+    block_length = L / real(sum(spacings))
+    h = block_length / real(block_size)
+    
+    idx = 1
+    x(idx) = 0.0
+    do b = 1, nblocks
+       if (periodic .and. (b == nblocks)) then
+          block_size = block_size - 1
+       end if
+       do i = 1, block_size
+          idx = idx + 1
+          x(idx) = x(idx - 1) + real(spacings(b)) * h
+       end do
+    end do
+    if (idx /= n) then
+       print *, "Error: didn't fill all interior nodes!"
+       error stop
+    end if
+
+    ! Add ghost points
+    if (.not. periodic) then
+       do i = 1, lghost
+          idx = 1 - i
+          x(idx) = x(idx + 1) - real(spacings(1)) * h
+       end do
+       do i = 1, rghost
+          idx = n + i
+          x(idx) = x(idx - 1) + real(spacings(nblocks)) * h
+       end do
+    else
+       idx = n + 1
+       x(idx) = L
+       do i = 2, rghost
+          idx = n + i
+          x(idx) = x(idx - 1) + (x(i) - x(i - 1))
+       end do
+
+       do i = 1, lghost
+          idx = 1 - i
+          x(idx) = x(idx + 1) - (x((n + 1) - (i - 1)) - (x((n + 1) - (i - 1) - 1)))
+       end do
+    end if
+
+    ! Self-check
+    if (any(x(1:n) > L) .or. any(x(1:n) < 0.0)) then
+       print *, "Error: grid exceeds specified length!"
+       print *, x
+       error stop
+    end if
+    if (abs(x(1)) > 0.0) then
+       print *, "Error: start point is wrong!"
+       print *, x(1)
+       error stop
+    end if
+
+    if (.not. periodic) then
+       idx = n
+    else
+       idx = n + 1
+    end if
+    if (abs(x(idx) - L) > (2 * epsilon(L) * L)) then
+       print *, "Error: end point is wrong!"
+       print *, "- Computed: ", x(idx)
+       print *, "- Expected: ", L
+       print *, "- Relative error: ", abs(x(idx) - L) / L
+       error stop
+    end if
+    
+  end subroutine build_grid
+  
+end program compute_derivatives
diff --git a/src/coeffs/nucfd_coeffs_a.f90 b/src/coeffs/nucfd_coeffs_a.f90
index 0d6ce03..98780ef 100644
--- a/src/coeffs/nucfd_coeffs_a.f90
+++ b/src/coeffs/nucfd_coeffs_a.f90
@@ -102,10 +102,10 @@ module subroutine coeff_a_components(h, numerator, numerator_corr, denominator,
     end select
 
     associate(beta => alpha) ! To match Gamet et al. (1999)
-      numerator = coeff_numerator(hm1, h0, hp2, hp2, beta)
+      numerator = coeff_numerator(hm1, h0, hp1, hp2, beta)
       numerator_corr = coeff_numerator_corr(hm1, h0, hp1, hp2, alpha, beta)
-      denominator = coeff_denominator(h0, hp1, hp2)
-      divisor = coeff_divisor(hm1, h0, hp1)
+      denominator = coeff_denominator(hm1, h0, hp1, hp2)
+      divisor = coeff_divisor(h0, hp1)
     end associate
 
   end subroutine coeff_a_components
@@ -145,29 +145,30 @@ pure real function coeff_numerator_corr(hm1, h0, hp1, hp2, alpha, beta)
   end function coeff_numerator_corr
 
   pure real function coeff_denominator(h0, hp1, hp2)
+  pure real function coeff_denominator(hm1, h0, hp1, hp2)
     !! Computes the denominator of the coefficient acting on f_{i+1}.
     !!
     !! Reduces to 3 h^3 when h=const. Dividing the numerator by this term yields the coefficient
     !! 14/9 when h=const, alpha=beta=1/3.
 
+    real, intent(in) :: hm1 
     real, intent(in) :: h0 
     real, intent(in) :: hp1 
     real, intent(in) :: hp2 
 
-    coeff_denominator = (3.0 * hp1 * ((h0 + hp1) / 2.0) * hp2)
+    coeff_denominator = (hp1 * (hm1 + h0 + hp1) * hp2)
   end function coeff_denominator
 
-  pure real function coeff_divisor(hm1, h0, hp1)
+  pure real function coeff_divisor(h0, hp1)
     !! Computes the non-uniform equivalent to 2h divisor of the coefficient acting on f_{i+1}.
     !!
     !! Dividing the coefficient by this term should reduce to (14/9)/(2h) when h=const,
     !! alpha=beta=1/3.
     
-    real, intent(in) :: hm1 
     real, intent(in) :: h0 
     real, intent(in) :: hp1 
 
-    coeff_divisor = (2.0 * ((hm1 + h0 + hp1) / 3.0))
+    coeff_divisor = (h0 + hp1)
   end function coeff_divisor
   
 end submodule nucfd_coeffs_a
diff --git a/src/coeffs/nucfd_coeffs_b.f90 b/src/coeffs/nucfd_coeffs_b.f90
index 2e104cc..c7b092b 100644
--- a/src/coeffs/nucfd_coeffs_b.f90
+++ b/src/coeffs/nucfd_coeffs_b.f90
@@ -103,8 +103,8 @@ module subroutine coeff_b_components(h, numerator, numerator_corr, denominator,
     associate(beta => alpha) ! To match Gamet et al. (1999)
       numerator = coeff_numerator(hm1, h0, hp1, hp2, alpha)
       numerator_corr = coeff_numerator_corr(hm1, h0, hp1, hp2, alpha, beta)
-      denominator = coeff_denominator(hm1, h0, hp1)
-      divisor = coeff_divisor(h0, hp1, hp2)
+      denominator = coeff_denominator(hm1, h0, hp1, hp2)
+      divisor = coeff_divisor(h0, hp1)
     end associate
 
   end subroutine coeff_b_components
@@ -143,7 +143,7 @@ pure real function coeff_numerator_corr(hm1, h0, hp1, hp2, alpha, beta)
                                                                     ! alpha = beta
   end function coeff_numerator_corr
 
-  pure real function coeff_denominator(hm1, h0, hp1)
+  pure real function coeff_denominator(hm1, h0, hp1, hp2)
     !! Computes the denominator of the coefficient acting on f_{i-1}.
     !!
     !! Reduces to 3 h^3 when h=const. Dividing the numerator by this term yields the coefficient
@@ -152,11 +152,12 @@ pure real function coeff_denominator(hm1, h0, hp1)
     real, intent(in) :: hm1
     real, intent(in) :: h0 
     real, intent(in) :: hp1
+    real, intent(in) :: hp2
     
-    coeff_denominator = 3.0 * hm1 * h0 * ((h0 + hp1) / 2.0)
+    coeff_denominator = hm1 * h0 * (h0 + hp1 + hp2)
   end function coeff_denominator
 
-  pure real function coeff_divisor(h0, hp1, hp2)
+  pure real function coeff_divisor(h0, hp1)
     !! Computes the non-uniform equivalent to 2h divisor of the coefficient acting on f_{i-1}.
     !!
     !! Dividing the coefficient by this term should reduce to (14/9)/(2h) when h=const,
@@ -164,9 +165,8 @@ pure real function coeff_divisor(h0, hp1, hp2)
 
     real, intent(in) :: h0
     real, intent(in) :: hp1
-    real, intent(in) :: hp2
 
-    coeff_divisor = 2.0 * (h0 + hp1 + hp2) / 3.0
+    coeff_divisor = (h0 + hp1) 
   end function coeff_divisor
   
 end submodule nucfd_coeffs_b
diff --git a/src/coeffs/nucfd_coeffs_c.f90 b/src/coeffs/nucfd_coeffs_c.f90
index 256a716..e90416a 100644
--- a/src/coeffs/nucfd_coeffs_c.f90
+++ b/src/coeffs/nucfd_coeffs_c.f90
@@ -106,8 +106,8 @@ module subroutine coeff_c_components(h, numerator, denominator, divisor)
 
     associate(beta => alpha) ! To match Gamet et al. (1999)
       numerator = coeff_numerator(hm1, h0, hp1, alpha, beta)
-      denominator = coeff_denominator(hm1, h0, hp1, hp2)
-      divisor = coeff_divisor(h0, hp1, hp2)
+      denominator = coeff_denominator(h0, hp1, hp2)
+      divisor = coeff_divisor(hm1, h0, hp1, hp2)
     end associate
     
   end subroutine coeff_c_components
@@ -129,33 +129,33 @@ pure real function coeff_numerator(hm1, h0, hp1, alpha, beta)
     
   end function coeff_numerator
   
-  pure real function coeff_denominator(hm1, h0, hp1, hp2)
+  pure real function coeff_denominator(h0, hp1, hp2)
     !! Computes the denominator of the coefficient acting on f_{i+2}.
     !!
     !! Reduces to 6 h^3 when h=const. Dividing the numerator by this term yields the coefficient
     !! 1/9 when h=const, alpha=beta=1/3.
 
-    real, intent(in) :: hm1
     real, intent(in) :: h0
     real, intent(in) :: hp1
     real, intent(in) :: hp2
 
-    coeff_denominator = 3.0 * hp2 * (hp1 + hp2) &
-         * ((hm1 + h0 + hp1 + hp2) / 4.0)
+    coeff_denominator = hp2 * (hp1 + hp2) &
+         * (h0 + hp1 + hp2)
     
   end function coeff_denominator
 
-  pure real function coeff_divisor(h0, hp1, hp2)
+  pure real function coeff_divisor(hm1, h0, hp1, hp2)
     !! Computes the non-uniform equivalent to 4h divisor of the coefficient acting on f_{i+2}.
     !!
     !! Dividing the coefficient by this term should reduce to (1/9)/(4h) when h=const,
     !! alpha=beta=1/3.
 
+    real, intent(in) :: hm1
     real, intent(in) :: h0
     real, intent(in) :: hp1
     real, intent(in) :: hp2
 
-    coeff_divisor = 4.0 * ((h0 + hp1 + hp2) / 3.0) 
+    coeff_divisor = (hm1 + h0 + hp1 + hp2) 
   end function coeff_divisor
   
 end submodule nucfd_coeffs_c
diff --git a/src/coeffs/nucfd_coeffs_d.f90 b/src/coeffs/nucfd_coeffs_d.f90
index 9dd0f80..b2d347d 100644
--- a/src/coeffs/nucfd_coeffs_d.f90
+++ b/src/coeffs/nucfd_coeffs_d.f90
@@ -106,8 +106,8 @@ module subroutine coeff_d_components(h, numerator, denominator, divisor)
 
     associate(beta => alpha) ! To match Gamet et al. (1999)
       numerator = coeff_numerator(h0, hp1, hp2, alpha, beta)
-      denominator = coeff_denominator(hm1, h0, hp1, hp2)
-      divisor = coeff_divisor(hm1, h0, hp1)
+      denominator = coeff_denominator(hm1, h0, hp1)
+      divisor = coeff_divisor(hm1, h0, hp1, hp2)
     end associate
     
   end subroutine coeff_d_components
@@ -129,7 +129,7 @@ pure real function coeff_numerator(h0, hp1, hp2, alpha, beta)
     
   end function coeff_numerator
   
-  pure real function coeff_denominator(hm1, h0, hp1, hp2)
+  pure real function coeff_denominator(hm1, h0, hp1)
     !! Computes the denominator of the coefficient acting on f_{f-2}.
     !!
     !! Reduces to 6 h^3 when h=const. Dividing the numerator by this term yields the coefficient
@@ -138,14 +138,13 @@ pure real function coeff_denominator(hm1, h0, hp1, hp2)
     real, intent(in) :: hm1
     real, intent(in) :: h0
     real, intent(in) :: hp1
-    real, intent(in) :: hp2
 
-    coeff_denominator = 3.0 * hm1 * (hp1 + hp2) &
-         * ((hm1 + h0 + hp1 + hp2) / 4.0)
+    coeff_denominator = hm1 * (hm1 + h0) &
+         * (hm1 + h0 + hp1)
     
   end function coeff_denominator
 
-  pure real function coeff_divisor(hm1, h0, hp1)
+  pure real function coeff_divisor(hm1, h0, hp1, hp2)
     !! Computes the non-uniform equivalent to 4h divisor of the coefficient acting on f_{i-2}.
     !!
     !! Dividing the coefficient by this term should reduce to -(1/9)/(4h) when h=const,
@@ -154,8 +153,9 @@ pure real function coeff_divisor(hm1, h0, hp1)
     real, intent(in) :: hm1
     real, intent(in) :: h0
     real, intent(in) :: hp1
+    real, intent(in) :: hp2
 
-    coeff_divisor = 4.0 * ((hm1 + h0 + hp1) / 3.0) 
+    coeff_divisor = (hm1 + h0 + hp1 + hp2)
   end function coeff_divisor
 
 end submodule nucfd_coeffs_d
diff --git a/src/nucfd_deriv_mod.f90 b/src/nucfd_deriv_mod.f90
index f9a946b..782aee9 100644
--- a/src/nucfd_deriv_mod.f90
+++ b/src/nucfd_deriv_mod.f90
@@ -28,6 +28,10 @@ subroutine deriv_rhs(f, stencil, x, dfdx)
     type(nucfd_stencil_points) :: stencil_coordinates
 
     real :: a, b, c, d, e
+
+    integer :: offset
+
+    offset = lbound(x, 1) - (-1)
     
     call create_stencil(5, 3, stencil_coordinates)
 
@@ -35,11 +39,7 @@ subroutine deriv_rhs(f, stencil, x, dfdx)
     type is(integer)
        select type(points => stencil_coordinates%stencil)
        type is(real)
-          points(-2) = x(indices(-2))
-          points(-1) = x(indices(-1))
-          points(+0) = x(indices(+0))
-          points(+1) = x(indices(+1))
-          points(+2) = x(indices(+2))
+          points(:) = x(indices(:) + offset)
        class default
           print *, "Error: Coordinate stencil is misallocated!"
           error stop
@@ -51,9 +51,9 @@ subroutine deriv_rhs(f, stencil, x, dfdx)
        d = coeff_d(stencil_coordinates)
        e = coeff_e(stencil_coordinates)
 
-       dfdx = a * (f(indices(+1)) + (b / a) * f(indices(-1))) &
-            + c * (f(indices(+2)) + (d / c) * f(indices(-2))) &
-            + e * f(indices(0))
+       dfdx = a * (f(indices(+1) + offset) + (b / a) * f(indices(-1) + offset)) &
+            + c * (f(indices(+2) + offset) + (d / c) * f(indices(-2) + offset)) &
+            + e * f(indices(0) + offset)
     class default
        print *, "Error: Index stencil is misallocated!"
        error stop
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index d2342c1..7467bdb 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -15,13 +15,14 @@ target_link_libraries(nucfd_test_flags INTERFACE
   nucfd_debug_flags
   nucfd_strict_flags)
 
-add_library(nucfd_test_framework nucfd_tests_mod.f90)
+add_library(nucfd_test_framework OBJECT
+  nucfd_tests_mod.f90)
 target_link_libraries(nucfd_test_framework nucfd_test_flags)
 
 function(define_test suite test_name)
   add_executable(${test_name} ${test_name}.f90)
   target_link_libraries(${test_name} nucfd)
-  target_link_libraries(${test_name} nucfd_test_framework)
+  target_link_libraries(${test_name} $<TARGET_OBJECTS:nucfd_test_framework>)
   target_link_libraries(${test_name} nucfd_test_flags)
   add_test(NAME ${suite}:${test_name} COMMAND ${test_name})
 endfunction()
diff --git a/tests/fd-schemes/differentiation_rules/CMakeLists.txt b/tests/fd-schemes/differentiation_rules/CMakeLists.txt
index 3ffc88a..d3234ea 100644
--- a/tests/fd-schemes/differentiation_rules/CMakeLists.txt
+++ b/tests/fd-schemes/differentiation_rules/CMakeLists.txt
@@ -9,6 +9,15 @@
 # SPDX-License-Identifier: BSD-3-Clause
 #
 
+add_library(diff_rules_utils OBJECT diff_rules_utils.f90)
+
 define_test(fd-schemes constant_function)
+add_test_libs(fd-schemes constant_function diff_rules_utils)
 define_test(fd-schemes linear_rising_function)
+add_test_libs(fd-schemes linear_rising_function diff_rules_utils)
 define_test(fd-schemes linear_falling_function)
+add_test_libs(fd-schemes linear_falling_function diff_rules_utils)
+define_test(fd-schemes quadratic_function)
+add_test_libs(fd-schemes quadratic_function diff_rules_utils)
+define_test(fd-schemes varmesh_constgrad)
+add_test_libs(fd-schemes varmesh_constgrad diff_rules_utils)
diff --git a/tests/fd-schemes/differentiation_rules/constant_function.f90 b/tests/fd-schemes/differentiation_rules/constant_function.f90
index 06ea223..5b19e77 100644
--- a/tests/fd-schemes/differentiation_rules/constant_function.f90
+++ b/tests/fd-schemes/differentiation_rules/constant_function.f90
@@ -11,50 +11,30 @@ program constant_function
   use nucfd_deriv
   
   use nucfd_tests
+  use diff_rules_utils
   
   implicit none
 
   integer :: n ! Grid size
   real :: L    ! Domain size
-  real :: h    ! Grid spacing
 
   real, dimension(:), allocatable :: x ! Grid
   real, dimension(:), allocatable :: f ! Function
   real :: dfdx ! Derivative
   real :: dgdx ! Derivative
 
-  ! Stencil
   type(nucfd_index_stencil) :: stencil
-  integer :: i
-  integer, parameter :: width = 5
-  integer, parameter :: centre = 3
+  type(test_setup) :: ts
   
   call initialise_suite("Constant function differentiation rules")
 
   n = 33
   L = 1.0
-  h = L / real(n - 1)
 
-  allocate(x(n))
-  allocate(f(n))
-  do i = 1, n
-    x(i) = real(i - 1) * h
-  end do
-
-  print *, "+++ Initialising stencil +++"
-  call create_stencil(width, centre, stencil)
-  i = 33 / 2
-  select type(indices => stencil%stencil)
-  type is(integer)
-     indices(-2) = i - 2
-     indices(-1) = i - 1
-     indices(+0) = i + 0
-     indices(+1) = i + 1
-     indices(+2) = i + 2
-  class default
-     print *, "Error: Index stencil is misallocated!"
-     error stop
-  end select
+  ts = initialise_test(n, L)
+  f = allocate_test_array(ts)
+  x = create_mesh_array(ts)
+  stencil = build_stencil(ts)
 
   f(:) = 1.0
   call deriv_rhs(f, stencil, x, dfdx)
diff --git a/tests/fd-schemes/differentiation_rules/diff_rules_utils.f90 b/tests/fd-schemes/differentiation_rules/diff_rules_utils.f90
new file mode 100644
index 0000000..94be47e
--- /dev/null
+++ b/tests/fd-schemes/differentiation_rules/diff_rules_utils.f90
@@ -0,0 +1,97 @@
+! tests/fd-schemes/differentiation_rules/diff_rules_utils.f90
+!
+!! Part of the fd-schemes test suite.
+!
+! SPDX-License-Identifier: BSD-3-Clause
+
+module diff_rules_utils
+  !! Utility module for differentiation rules test subsuite.
+
+  use nucfd_types
+  
+  implicit none
+
+  type test_setup
+     integer :: n
+     real :: L
+     real :: h
+  end type test_setup
+  
+  private
+  public :: test_setup
+  public :: initialise_test
+  public :: allocate_test_array
+  public :: create_mesh_array
+  public :: build_stencil
+  
+  integer, parameter :: width = 5
+  integer, parameter :: centre = 3
+  
+contains
+
+  type(test_setup) function initialise_test(n, L)
+
+    integer, intent(in) :: n
+    real, intent(in) :: L
+
+    real :: h
+
+    h = L / real(n - 1)
+    initialise_test = test_setup(n, L, h)
+
+  end function initialise_test
+
+  function allocate_test_array(ts) result(arr)
+
+    type(test_setup), intent(in) :: ts
+    real, dimension(:), allocatable :: arr
+
+    allocate(arr(ts%n))
+    
+  end function allocate_test_array
+
+  function create_mesh_array(ts) result(x)
+
+    type(test_setup), intent(in) :: ts
+    real, dimension(:), allocatable :: x
+
+    integer :: i
+    
+    x = allocate_test_array(ts)
+
+    associate(n => ts%n, h => ts%h)
+      do i = 1, n
+         if (i < n / 2) then
+            x(i) = real(i - 1) * h
+         else
+            x(i) = x(i - 1) + h / 2.0
+         end if
+      end do
+    end associate
+
+  end function create_mesh_array
+
+  type(nucfd_index_stencil) function build_stencil(ts)
+
+    type(test_setup), intent(in) :: ts
+
+    integer :: i
+    
+    print *, "+++ Initialising stencil +++"
+    call create_stencil(width, centre, build_stencil)
+    i = ts%n / 2
+    select type(indices => build_stencil%stencil)
+    type is(integer)
+       indices(-2) = i - 2
+       indices(-1) = i - 1
+       indices(+0) = i + 0
+       indices(+1) = i + 1
+       indices(+2) = i + 2
+    class default
+       print *, "Error: Index stencil is misallocated!"
+       error stop
+    end select
+
+  end function build_stencil
+  
+end module diff_rules_utils
diff --git a/tests/fd-schemes/differentiation_rules/linear_falling_function.f90 b/tests/fd-schemes/differentiation_rules/linear_falling_function.f90
index eab0827..fd16fc8 100644
--- a/tests/fd-schemes/differentiation_rules/linear_falling_function.f90
+++ b/tests/fd-schemes/differentiation_rules/linear_falling_function.f90
@@ -11,55 +11,37 @@ program linear_falling_function
   use nucfd_deriv
   
   use nucfd_tests
+  use diff_rules_utils
   
   implicit none
 
   integer :: n ! Grid size
   real :: L    ! Domain size
-  real :: h    ! Grid spacing
 
   real, dimension(:), allocatable :: x ! Grid
   real, dimension(:), allocatable :: f ! Function
   real :: dfdx ! Derivative
   real :: dgdx ! Derivative
 
-  ! Stencil
   type(nucfd_index_stencil) :: stencil
+  type(test_setup) :: ts
+
   integer :: i
-  integer, parameter :: width = 5
-  integer, parameter :: centre = 3
   
   call initialise_suite("Linearly decreasing function differentiation rules")
 
   n = 33
   L = 1.0
-  h = L / real(n - 1)
-
-  allocate(x(n))
-  allocate(f(n))
-  do i = 1, n
-    x(i) = real(i - 1) * h
-  end do
 
-  print *, "+++ Initialising stencil +++"
-  call create_stencil(width, centre, stencil)
-  i = 33 / 2
-  select type(indices => stencil%stencil)
-  type is(integer)
-     indices(-2) = i - 2
-     indices(-1) = i - 1
-     indices(+0) = i + 0
-     indices(+1) = i + 1
-     indices(+2) = i + 2
-  class default
-     print *, "Error: Index stencil is misallocated!"
-     error stop
-  end select
+  ts = initialise_test(n, L)
+  f = allocate_test_array(ts)
+  x = create_mesh_array(ts)
+  stencil = build_stencil(ts)
 
   f(1) = 0.0
   dfdx = 1.0
   do i = 2, n
-     f(i) = f(i - 1) - h * dfdx
+     f(i) = f(i - 1) - ts%h * dfdx
   end do
 
   call deriv_rhs(f, stencil, x, dfdx)
diff --git a/tests/fd-schemes/differentiation_rules/linear_rising_function.f90 b/tests/fd-schemes/differentiation_rules/linear_rising_function.f90
index 4400537..05d28b2 100644
--- a/tests/fd-schemes/differentiation_rules/linear_rising_function.f90
+++ b/tests/fd-schemes/differentiation_rules/linear_rising_function.f90
@@ -11,55 +11,37 @@ program linear_rising_function
   use nucfd_deriv
   
   use nucfd_tests
+  use diff_rules_utils
   
   implicit none
 
   integer :: n ! Grid size
   real :: L    ! Domain size
-  real :: h    ! Grid spacing
 
   real, dimension(:), allocatable :: x ! Grid
   real, dimension(:), allocatable :: f ! Function
   real :: dfdx ! Derivative
   real :: dgdx ! Derivative
 
-  ! Stencil
   type(nucfd_index_stencil) :: stencil
+  type(test_setup) :: ts
+
   integer :: i
-  integer, parameter :: width = 5
-  integer, parameter :: centre = 3
   
   call initialise_suite("Linearly increasing function differentiation rules")
 
   n = 33
   L = 1.0
-  h = L / real(n - 1)
-
-  allocate(x(n))
-  allocate(f(n))
-  do i = 1, n
-    x(i) = real(i - 1) * h
-  end do
 
-  print *, "+++ Initialising stencil +++"
-  call create_stencil(width, centre, stencil)
-  i = 33 / 2
-  select type(indices => stencil%stencil)
-  type is(integer)
-     indices(-2) = i - 2
-     indices(-1) = i - 1
-     indices(+0) = i + 0
-     indices(+1) = i + 1
-     indices(+2) = i + 2
-  class default
-     print *, "Error: Index stencil is misallocated!"
-     error stop
-  end select
+  ts = initialise_test(n, L)
+  f = allocate_test_array(ts)
+  x = create_mesh_array(ts)
+  stencil = build_stencil(ts)
 
   f(1) = 0.0
   dfdx = 1.0
   do i = 2, n
-     f(i) = f(i - 1) + h * dfdx
+     f(i) = f(i - 1) + ts%h * dfdx
   end do
   call deriv_rhs(f, stencil, x, dfdx)
   call deriv_rhs(f + 1.0, stencil, x, dgdx)
diff --git a/tests/fd-schemes/differentiation_rules/quadratic_function.f90 b/tests/fd-schemes/differentiation_rules/quadratic_function.f90
new file mode 100644
index 0000000..218be94
--- /dev/null
+++ b/tests/fd-schemes/differentiation_rules/quadratic_function.f90
@@ -0,0 +1,61 @@
+! tests/fd-schemes/differentiation_rules/quadratic_function.f90
+!
+!! Part of the fd-schemes test suite.
+!
+! SPDX-License-Identifier: BSD-3-Clause
+
+program quadratic_function
+  !! Tests that rules of differentiation are respected for quadratic functions.
+
+  use nucfd_types
+  use nucfd_deriv
+  
+  use nucfd_tests
+  use diff_rules_utils
+  
+  implicit none
+
+  integer :: n ! Grid size
+  real :: L    ! Domain size
+
+  real, dimension(:), allocatable :: x ! Grid
+  real, dimension(:), allocatable :: f ! Function
+  real :: dfdx ! Derivative
+  real :: dgdx ! Derivative
+
+  type(nucfd_index_stencil) :: stencil
+  type(test_setup) :: ts
+
+  integer :: i
+  
+  call initialise_suite("Quadratic function differentiation rules")
+
+  n = 33
+  L = 1.0
+
+  ts = initialise_test(n, L)
+  f = allocate_test_array(ts)
+  x = create_mesh_array(ts)
+  stencil = build_stencil(ts)
+
+  do i = 1, n
+     f(i) = x(i)**2
+  end do
+  call deriv_rhs(f, stencil, x, dfdx)
+  call deriv_rhs(f + 1.0, stencil, x, dgdx)
+  call test_report("Shifted derivative +, linear+ f", check_scalar(dgdx, dfdx))
+  call deriv_rhs(f - 1.0, stencil, x, dgdx)
+  call test_report("Shifted derivative -, linear+ f", check_scalar(dgdx, dfdx))
+  call deriv_rhs(f, stencil, x, dgdx)
+  call test_report("Shifted derivative 0, linear+ f", check_scalar(dgdx, dfdx))
+  call deriv_rhs(2.0 * f, stencil, x, dgdx)
+  call test_report("Scaled derivative 2x, linear+ f", check_scalar(dgdx, 2.0 * dfdx))
+  call deriv_rhs(-f, stencil, x, dgdx)
+  call test_report("Scaled derivative -1, linear+ f", check_scalar(dgdx, -dfdx))
+  
+  deallocate(x)
+  deallocate(f)
+  
+  call finalise_suite()
+  
+end program quadratic_function
diff --git a/tests/fd-schemes/differentiation_rules/varmesh_constgrad.f90 b/tests/fd-schemes/differentiation_rules/varmesh_constgrad.f90
new file mode 100644
index 0000000..aea4db2
--- /dev/null
+++ b/tests/fd-schemes/differentiation_rules/varmesh_constgrad.f90
@@ -0,0 +1,87 @@
+! tests/fd-schemes/differentiation_rules/varmesh_constgrad.f90
+!
+!! Part of the fd-schemes test suite.
+!
+! SPDX-License-Identifier: BSD-3-Clause
+
+program linear_rising_function
+  !! Tests that the derivative of a linear function remains constant on a mesh with local
+  !! refinement.
+
+  use nucfd_types
+  use nucfd_deriv
+  
+  use nucfd_tests
+  use diff_rules_utils
+  
+  implicit none
+
+  integer :: n ! Grid size
+  real :: L    ! Domain size
+
+  real, dimension(:), allocatable :: x ! Grid
+  real, dimension(:), allocatable :: f ! Function
+  real :: dfdx ! Derivative
+  real :: dgdx ! Derivative
+
+  type(nucfd_index_stencil) :: stencil
+  type(test_setup) :: ts
+
+  integer :: i, idx
+
+  logical :: passing
+  
+  call initialise_suite("Constant gradient on locally refined mesh")
+
+  n = 33
+  L = 1.0
+
+  ts = initialise_test(n, L)
+  f = allocate_test_array(ts)
+  x = create_mesh_array(ts)
+  stencil = build_stencil(ts)
+
+  f(1) = 0.0
+  dfdx = 1.0
+  do i = 2, n
+     f(i) = f(i - 1) + (x(i) - x(i - 1)) * dfdx
+  end do
+
+  i = 1
+  call centre_stencil(i, stencil)
+  call deriv_rhs(f, stencil, x, dfdx)
+
+  passing = .true.
+  do i = 2, n - 4
+     idx = i + 2
+     print *, idx, x(idx - 1) - x(idx - 2), x(idx) - x(idx - 1), x(idx + 1) - x(idx), x(idx + 2) - x(idx + 1)
+     call centre_stencil(i, stencil)
+     call deriv_rhs(f, stencil, x, dgdx)
+     passing = passing .and. check_scalar(dgdx, dfdx)
+  end do
+  call test_report("Gradient remains constant", passing)
+  
+  deallocate(x)
+  deallocate(f)
+  
+  call finalise_suite()
+
+contains
+
+  subroutine centre_stencil(idx, stencil)
+    !! Move a stencil to be centred at specified index.
+
+    integer, intent(in) :: idx                          !! The index
+    type(nucfd_index_stencil), intent(inout) :: stencil !! The stencil to be moved
+
+    integer :: shift
+
+    select type(indices => stencil%stencil)
+    type is(integer)
+       shift = idx - indices(0)
+       indices(:) = indices(:) + shift
+    end select
+    
+  end subroutine centre_stencil
+ 
+end program linear_rising_function
diff --git a/tests/fd-schemes/verify_coeffs/CMakeLists.txt b/tests/fd-schemes/verify_coeffs/CMakeLists.txt
index f95f68c..26e3d3f 100644
--- a/tests/fd-schemes/verify_coeffs/CMakeLists.txt
+++ b/tests/fd-schemes/verify_coeffs/CMakeLists.txt
@@ -14,3 +14,4 @@ define_test(fd-schemes verify_coeff_b)
 define_test(fd-schemes verify_coeff_c)
 define_test(fd-schemes verify_coeff_d)
 define_test(fd-schemes verify_coeff_e)
+define_test(fd-schemes verify_coeff_pairsums)
diff --git a/tests/fd-schemes/verify_coeffs/verify_coeff_a.f90 b/tests/fd-schemes/verify_coeffs/verify_coeff_a.f90
index 9c3e5b4..f12f589 100644
--- a/tests/fd-schemes/verify_coeffs/verify_coeff_a.f90
+++ b/tests/fd-schemes/verify_coeffs/verify_coeff_a.f90
@@ -4,7 +4,7 @@
 !
 ! SPDX-License-Identifier: BSD-3-Clause
 
-program verify_coeff_b
+program verify_coeff_a
   !! Tests the computation of finite difference coefficient for non-uniform grids acting on f_{i+1}.
 
   use nucfd_types
@@ -21,13 +21,10 @@ program verify_coeff_b
   type(nucfd_stencil_points) :: stencil ! Stencil of grid points.
   
   real :: a
-  real, parameter :: aref = 14.0 / 9.0
 
-  real :: numerator, numerator_corr, denominator, divisor
-  real, parameter :: numerator_f1ref = 14.0 / 3.0
-  real, parameter :: numerator_corr_f1ref = 0.0
-  real, parameter :: denominator_f1ref = 3.0
-  real, parameter :: divisor_f1ref = 2.0
+  integer :: scale
+  real :: aref
+  real :: numerator_f1ref, numerator_corr_f1ref, denominator_f1ref, divisor_f1ref
   
   call initialise_suite("Verify coefficient A")
 
@@ -36,26 +33,90 @@ program verify_coeff_b
   h = L / real(n - 1)
 
   call create_stencil(5, 3, stencil)
-  select type(points => stencil%stencil)
-  type is(real)
-     points(:) = 0.0
-     points(-2) = -2.0 * h
-     points(-1) = -1.0 * h
-     points(0) = 0.0
-     points(1) = +1.0 * h
-     points(2) = +2.0 * h
-  class default
-     print *, "Error: Coordinate stencil is misallocated!"
-  end select
-
-  call coeff_a_components(points_to_deltas(stencil), numerator, numerator_corr, denominator, divisor)
-  call test_report("Coefficient A numerator", check_scalar(numerator, numerator_f1ref * (h**3)))
-  call test_report("Coefficient A numerator correction", check_scalar(numerator_corr, numerator_corr_f1ref * (h**3)))
-  call test_report("Coefficient A denominator", check_scalar(denominator, denominator_f1ref * (h**3)))
-  call test_report("Coefficient A divisor", check_scalar(divisor, divisor_f1ref * h))
-  a = coeff_a(stencil)
-  call test_report("Coefficient A", check_scalar(a, aref / (2.0 * h)))
 
+  scale = 1
+  call setup_stencil(scale, scale, scale, scale)
+  call compute_uniform_ref_coeffs(scale)
+  call test_coefficient("Uniform", stencil, aref, numerator_f1ref, numerator_corr_f1ref, &
+       denominator_f1ref, divisor_f1ref)
+
+  !! Check coefficient is independent of scale
+  scale = 2
+  call setup_stencil(scale, scale, scale, scale)
+  call compute_uniform_ref_coeffs(scale)
+  call test_coefficient("Uniform (2h)", stencil, aref, numerator_f1ref, numerator_corr_f1ref, &
+       denominator_f1ref, divisor_f1ref)
+  
   call finalise_suite()
 
-end program verify_coeff_b
+contains
+
+  subroutine setup_stencil(sm2, sm1, sp1, sp2)
+
+    integer, intent(in) :: sm2
+    integer, intent(in) :: sm1
+    integer, intent(in) :: sp1
+    integer, intent(in) :: sp2
+
+    select type(points => stencil%stencil)
+    type is(real)
+       points(:) = 0.0
+       points(-2) = -(sm1 + sm2) * h
+       points(-1) = -sm1 * h
+       points(0) = 0.0
+       points(1) = +sp1 * h
+       points(2) = +(sp1 + sp2) * h
+    class default
+       print *, "Error: Coordinate stencil is misallocated!"
+    end select
+    
+  end subroutine setup_stencil
+
+  subroutine compute_uniform_ref_coeffs(scale)
+
+    integer, intent(in) :: scale
+    
+    aref = 14.0 / 9.0 / (2.0 * (scale * h))
+    numerator_f1ref = 14.0 / 3.0 * ((scale * h)**3)
+    numerator_corr_f1ref = 0.0 * ((scale * h)**3)
+    denominator_f1ref = 3.0 * ((scale * h)**3)
+    divisor_f1ref = 2.0 * (scale * h)
+  end subroutine compute_uniform_ref_coeffs
+  
+  subroutine test_coefficient(grid_desc, stencil, aref, numerator_f1ref, numerator_corr_f1ref, &
+       denominator_f1ref, divisor_f1ref)
+
+    character(len=*), intent(in) :: grid_desc
+    type(nucfd_stencil_points), intent(in) :: stencil
+    real, intent(in) :: aref
+    real, intent(in) :: numerator_f1ref
+    real, intent(in) :: numerator_corr_f1ref
+    real, intent(in) :: denominator_f1ref
+    real, intent(in) :: divisor_f1ref
+
+    character(len=:), allocatable :: prefix
+    real :: numerator, numerator_corr, denominator, divisor
+
+    prefix = "("//grid_desc//") "
+    
+    call coeff_a_components(points_to_deltas(stencil), numerator, numerator_corr, &
+         denominator, divisor)
+    call test_report(prefix//"Coefficient A numerator", &
+         check_scalar(numerator, numerator_f1ref))
+    call test_report(prefix//"Coefficient A numerator correction", &
+         check_scalar(numerator_corr, numerator_corr_f1ref))
+    call test_report(prefix//"Coefficient A denominator", &
+         check_scalar(denominator, denominator_f1ref))
+    call test_report(prefix//"Coefficient A divisor", &
+         check_scalar(divisor, divisor_f1ref))
+    call test_report(prefix//"Coefficient A den * div", &
+         check_scalar(denominator * divisor, &
+         (denominator_f1ref * divisor_f1ref)))
+    
+    a = coeff_a(stencil)
+    call test_report("Coefficient A", &
+         check_scalar(a, aref))
+
+  end subroutine test_coefficient
+  
+end program verify_coeff_a
diff --git a/tests/fd-schemes/verify_coeffs/verify_coeff_e.f90 b/tests/fd-schemes/verify_coeffs/verify_coeff_e.f90
index 9f5f75a..55409e0 100644
--- a/tests/fd-schemes/verify_coeffs/verify_coeff_e.f90
+++ b/tests/fd-schemes/verify_coeffs/verify_coeff_e.f90
@@ -43,7 +43,7 @@ program verify_coeff_e
   end select
 
   e = coeff_e(stencil)
-  call test_report("Coefficient E", check_scalar(e, eref / (4.0 * h)))
+  call test_report("Coefficient E", check_scalar(e, eref))
 
   call finalise_suite()
 
diff --git a/tests/fd-schemes/verify_coeffs/verify_coeff_pairsums.f90 b/tests/fd-schemes/verify_coeffs/verify_coeff_pairsums.f90
new file mode 100644
index 0000000..609c39a
--- /dev/null
+++ b/tests/fd-schemes/verify_coeffs/verify_coeff_pairsums.f90
@@ -0,0 +1,55 @@
+! tests/fd-schemes/verify_coeffs/verify_coeff_pairsums.f90
+!
+!! Part of the fd-schemes test suite.
+!
+! SPDX-License-Identifier: BSD-3-Clause
+
+program verify_coeff_pairsums
+  !! Tests the computation of finite difference coefficients for the pairs f_{i+/-1} and f_{i+/-2}.
+
+  use nucfd_types
+  use nucfd_coeffs
+
+  use nucfd_tests
+  
+  implicit none
+  
+  real :: n ! Mesh size
+  real :: L ! Domain size
+  real :: h ! Average grid spacing.
+
+  type(nucfd_stencil_points) :: stencil ! Stencil of grid points.
+  
+  real :: a, b, c, d
+  real, parameter :: abref = 0.0, cdref = 0.0
+  
+  call initialise_suite("Verify coefficient pairsums")
+
+  n = 128
+  L = 1.0
+  h = L / real(n - 1)
+
+  call create_stencil(5, 3, stencil)
+  select type(points => stencil%stencil)
+  type is(real)
+     points(:) = 0.0
+     points(-2) = -2.0 * h
+     points(-1) = -1.0 * h
+     points(0) = 0.0
+     points(1) = +1.0 * h
+     points(2) = +2.0 * h
+  class default
+     print *, "Error: Coordinate stencil is misallocated!"
+  end select
+
+  a = coeff_a(stencil)
+  b = coeff_b(stencil)
+  call test_report("Coefficient A+B", check_scalar(a + b, abref))
+
+  c = coeff_c(stencil)
+  d = coeff_d(stencil)
+  call test_report("Coefficient A+B", check_scalar(c + d, cdref))
+
+  call finalise_suite()
+
+end program verify_coeff_pairsums
diff --git a/tests/tridsolver/CMakeLists.txt b/tests/tridsolver/CMakeLists.txt
index 56cb4b5..2a9de12 100644
--- a/tests/tridsolver/CMakeLists.txt
+++ b/tests/tridsolver/CMakeLists.txt
@@ -9,15 +9,16 @@
 # SPDX-License-Identifier: BSD-3-Clause
 #
 
-add_library(trid_test_utils tridsol_test_utils_mod.f90)
+add_library(trid_test_utils OBJECT
+  tridsol_test_utils_mod.f90)
 target_link_libraries(trid_test_utils nucfd_test_flags)
 
 define_test(tridsolver system_11_symm nucfd_test_framework)
 define_test(tridsolver system_11_anti-symm)
-add_test_libs(tridsolver system_11_symm trid_test_utils)
-add_test_libs(tridsolver system_11_anti-symm trid_test_utils)
+add_test_libs(tridsolver system_11_symm $<TARGET_OBJECTS:trid_test_utils>)
+add_test_libs(tridsolver system_11_anti-symm $<TARGET_OBJECTS:trid_test_utils>)
 
 define_test(tridsolver system_00_symm nucfd_test_framework)
 define_test(tridsolver system_00_anti-symm)
-add_test_libs(tridsolver system_00_symm trid_test_utils)
-add_test_libs(tridsolver system_00_anti-symm trid_test_utils)
+add_test_libs(tridsolver system_00_symm $<TARGET_OBJECTS:trid_test_utils>)
+add_test_libs(tridsolver system_00_anti-symm $<TARGET_OBJECTS:trid_test_utils>)