diff --git a/CHANGELOG.md b/CHANGELOG.md
index 2eff767..c180acf 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -6,12 +6,13 @@ This is the change log for `NuCFD`, it is based on the format suggested by [Keep
 
 ### Added
 
-- Added Cray compiler support [42b7070].
-- Added Intel compiler support [2d3ce07].
 - Added `FORD` documentation system [c780fe1].
 - Added tridiagonal solver with support for periodic problems [f5d254f].
 - Added simple library to support writing tests [a5814c8].
 - Added build and test system based on `cmake` [8d56f30].
+  - Added Cray compiler support [42b7070].
+  - Added Intel compiler support [2d3ce07].
+  - Set minimum `cmake` version `v3.25.2` to support Fortran submodules [e8c4aba].
 
 ### Changed
 ### Deprecated
diff --git a/CMakeLists.txt b/CMakeLists.txt
index a658090..2f62c41 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -2,16 +2,23 @@
 #
 ## Description
 #
+# Main CMakeLists file for the NuCFD project.
+#
 ## LICENSE
 #
 # SPDX-License-Identifier: BSD-3-Clause
 #
 
 ## Basic configuration
-cmake_minimum_required(VERSION 3.16)
+cmake_minimum_required(VERSION 3.25.2)
 project(NuCFD)
 enable_language(Fortran)
 
+option(BUILD_SHARED_LIBS "Enable building as a shared library" OFF)
+if(BUILD_SHARED_LIBS)
+  set(CMAKE_POSITION_INDEPENDENT_CODE true)
+endif()
+
 include(GNUInstallDirs)
 set(LIBRARY_OUTPUT_PATH ${PROJECT_BINARY_DIR}/${CMAKE_INSTALL_LIBDIR})
 set(EXECUTABLE_OUTPUT_PATH ${PROJECT_BINARY_DIR}/${CMAKE_INSTALL_BINDIR})
@@ -23,9 +30,9 @@ set(icc_like_fc "$<COMPILE_LANG_AND_ID:Fortran,Intel>")
 set(cray_like_fc "$<COMPILE_LANG_AND_ID:Fortran,Cray>")
 add_library(nucfd_compiler_flags INTERFACE)
 target_compile_options(nucfd_compiler_flags INTERFACE
-  "$<${gcc_like_fc}:-std=f2008;-Wall;-Wextra;-Wpedantic;-Warray-bounds;-Wimplicit-interface;-Wimplicit-procedure;-fimplicit-none>"
-  "$<${icc_like_fc}:-std08;-warn;all>"
-  "$<${cray_like_fc}:-en;-emf;-eI>")
+  "$<${gcc_like_fc}:-cpp;-std=f2018;-Wall;-Wextra;-Wpedantic;-Warray-bounds;-Wimplicit-interface;-Wimplicit-procedure;-fimplicit-none>"
+  "$<${icc_like_fc}:-fpp;-std18;-warn;all>"
+  "$<${cray_like_fc}:-cpp;-hfp0;-en;-emf;-eI>") # Anything greater than fp0 fails tests!
 add_library(nucfd_debug_flags INTERFACE)
 target_compile_options(nucfd_debug_flags INTERFACE
   "$<${gcc_like_fc}:-g;-Og;-ffpe-trap=invalid,zero;-fcheck=bounds;-fbacktrace>"
@@ -35,48 +42,16 @@ add_library(nucfd_strict_flags INTERFACE)
 target_compile_options(nucfd_strict_flags INTERFACE
   "$<${gcc_like_fc}:-Werror;-Wno-error=integer-division>"
   "$<${cray_like_fc}:-eN>")
+if (CMAKE_BUILD_TYPE MATCHES "Debug")
+  target_link_libraries(nucfd_compiler_flags INTERFACE nucfd_debug_flags)
+endif()
 
 ## Build nucfd library
-add_library(nucfd
-  src/nucfd_trid_solver_mod.f90
-  src/nucfd_coeffs_mod.f90)
-target_link_libraries(nucfd nucfd_compiler_flags)
+add_subdirectory(src)
 
 ## Testing
 enable_testing()
-
-add_library(nucfd_test_flags INTERFACE)
-target_link_libraries(nucfd_test_flags INTERFACE
-  nucfd_compiler_flags
-  nucfd_debug_flags
-  nucfd_strict_flags)
-
-add_library(nucfd_test_framework tests/nucfd_tests_mod.f90)
-target_link_libraries(nucfd_test_framework nucfd_test_flags)
-
-function(define_test suite test_name)
-  add_executable(${test_name} tests/${suite}/${test_name}.f90)
-  target_link_libraries(${test_name} nucfd)
-  target_link_libraries(${test_name} nucfd_test_framework)
-  target_link_libraries(${test_name} nucfd_test_flags)
-  add_test(NAME ${suite}:${test_name} COMMAND ${test_name})
-endfunction()
-
-function(add_test_libs suite test_name lib)
-  target_link_libraries(${test_name} ${lib})
-endfunction()
-
-# Tidiagonal solver test suite
-add_library(trid_test_utils tests/tridsolver/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)
-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_subdirectory(tests)
 
 ## Documentation
 add_custom_target(doc ford
diff --git a/README.md b/README.md
index 886b87c..b92451c 100644
--- a/README.md
+++ b/README.md
@@ -4,36 +4,49 @@ difference schemes on non-uniform grids.
 ## Building
 
 The `NuCFD` build system is generated by `cmake`, a default configuration can be created by running
-``
+```
 cmake -B build
-``
+```
 from the root directory.
 Once the build system has been generated running
-``
+```
 make -C build
-``
+```
 will build the library.
 The build can be configured using the `ccmake` tool, however currently only `gfortran` is supported.
 
+Note that this project uses Fortran2008 submodules, support for these in CMake requires CMake v3.25.2
+or above.
+If you already have an older CMake it can be relatively easily updated for the local user by cloning
+the repository and checking out tag `v3.25.2` (or later), CMake iteslf can then be built using a
+fairly standard process:
+```
+cmake -B build -DCMAKE_INSTALL_PREFIX=/path/to/install
+# Any configuration required
+make -C build && make -C build install
+export PATH=/path/to/install/bin:${PATH}
+```
+you should then be able to use your new CMake to configure the `NuCFD` build.
+
 ## Testing
 
 After building the library it can be tested using `ctest`.
 From the root directory run
-``
+```
 make -C build test
-``
+```
 to launch the tests.
 If any test fail more detailed output can be shown by running `ctest` verbosely
-``
+```
 cd build
 ctest --verbose
-``
+```
 which will print any output from the tests, add `--rerun-failed` to only repeat failing tests.
 
 ## Documentation
 
 Documentation can be generated using the FORD tool by running
-``
+```
 make -C build/ doc
-``
+```
 the generated documentation can be viewed in a web browser by opening `./manual/index.html`.
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
new file mode 100644
index 0000000..9efc08b
--- /dev/null
+++ b/src/CMakeLists.txt
@@ -0,0 +1,24 @@
+## src/CMakeLists.txt
+#
+## Description
+#
+# CMakeLists file for building the NuCFD library.
+#
+## LICENSE
+#
+# SPDX-License-Identifier: BSD-3-Clause
+#
+
+add_library(nucfd_types OBJECT
+  nucfd_types_mod.f90
+)
+target_link_libraries(nucfd_types PRIVATE nucfd_compiler_flags)
+
+add_subdirectory(coeffs)
+
+add_library(nucfd
+  nucfd_trid_solver_mod.f90
+  nucfd_deriv_mod.f90)
+target_link_libraries(nucfd PRIVATE $<TARGET_OBJECTS:nucfd_types>)
+target_link_libraries(nucfd PRIVATE $<TARGET_OBJECTS:nucfd_coeffs>)
+target_link_libraries(nucfd PUBLIC nucfd_compiler_flags)
diff --git a/src/coeffs/CMakeLists.txt b/src/coeffs/CMakeLists.txt
new file mode 100644
index 0000000..d0b4672
--- /dev/null
+++ b/src/coeffs/CMakeLists.txt
@@ -0,0 +1,21 @@
+## src/coeffs/CMakeLists.txt
+#
+## Description
+#
+# CMakeLists file for building the coefficients module of NuCFD.
+#
+## LICENSE
+#
+# SPDX-License-Identifier: BSD-3-Clause
+#
+
+add_library(nucfd_coeffs OBJECT
+  nucfd_coeffs_mod.f90
+  nucfd_coeffs_a.f90
+  nucfd_coeffs_b.f90
+  nucfd_coeffs_c.f90
+  nucfd_coeffs_d.f90
+  nucfd_coeffs_e.f90
+)
+target_link_libraries(nucfd_coeffs PRIVATE $<TARGET_OBJECTS:nucfd_types>)
+target_link_libraries(nucfd_coeffs PRIVATE nucfd_compiler_flags)
diff --git a/src/coeffs/nucfd_coeffs_a.f90 b/src/coeffs/nucfd_coeffs_a.f90
new file mode 100644
index 0000000..0d6ce03
--- /dev/null
+++ b/src/coeffs/nucfd_coeffs_a.f90
@@ -0,0 +1,173 @@
+! src/coeffs/nucfd_coeff_a.f90
+!
+!! Implements coefficient A of the nucfd_coeffs module.
+!
+! SPDX-License-Identifier: BSD-3-Clause
+
+submodule (nucfd_coeffs) nucfd_coeffs_a
+  !! Submodule defining the coefficient acting on f_{i+1} for compact finite difference schemes on
+  !! non-uniform grids.
+
+  implicit none
+
+contains
+
+  module real function coeff_a_points(x)
+    !! Compute the coefficient acting on f_{i+1} of the finite diference given a stencil of points.
+
+    type(nucfd_stencil_points), intent(in) :: x !! Stencil of points for the finite
+                                                !! difference.
+
+    type(nucfd_stencil_deltas) :: h  ! Stencil of grid spacings for the finite difference.
+
+#ifndef NDEBUG
+    if (size(x%stencil) /= 5) then
+       print *, "Error@coeff_a_points: expecting width 5 stencil, received width = ", &
+            size(x%stencil)
+       error stop
+    end if
+    if (1 - lbound(x%stencil, 1) /= 3) then
+       print *, "Error@coeff_a_points: expecting centre 3 stencil, received centre = ", &
+            1 - lbound(x%stencil, 1)
+       print *, size(x%stencil), lbound(x%stencil), ubound(x%stencil)
+       error stop
+    end if
+#endif
+    
+    h = points_to_deltas(x)
+    coeff_a_points = coeff_a(h)
+    
+  end function coeff_a_points
+
+  module real function coeff_a_deltas(h)
+    !! Compute the coefficient acting on f_{i+1} of the finite diference given a stencil of grid
+    !! spacings.
+
+    type(nucfd_stencil_deltas), intent(in) :: h !! Stencil of grid spacings for the finite difference.
+
+    real :: numerator, numerator_corr, denominator, divisor
+
+#ifndef NDEBUG
+    if (size(h%stencil) /= 4) then
+       print *, "Error@coeff_a_deltas: expecting width 4 stencil, received width = ", &
+            size(h%stencil)
+       error stop
+    end if
+    if (1 - lbound(h%stencil, 1) /= 2) then
+       print *, "Error@coeff_a_deltas: expecting centre 2 stencil, received centre = ", &
+            1 - lbound(h%stencil, 1)
+       print *, size(h%stencil), lbound(h%stencil), ubound(h%stencil)
+       error stop
+    end if
+#endif
+
+    call coeff_a_components(h, numerator, numerator_corr, denominator, divisor)
+    coeff_a_deltas = ((numerator + numerator_corr) / denominator) / divisor
+    
+  end function coeff_a_deltas
+
+  module subroutine coeff_a_components(h, numerator, numerator_corr, denominator, divisor)
+
+    type(nucfd_stencil_deltas), intent(in) :: h !! Stencil of grid spacings for the finite difference.
+
+    real, intent(out) :: numerator
+    real, intent(out) :: numerator_corr
+    real, intent(out) :: denominator
+    real, intent(out) :: divisor
+    
+    real :: hm1, h0, hp1, hp2 ! Grid deltas at i -2, -1, 0, +1, +2
+
+#ifndef NDEBUG
+    if (size(h%stencil) /= 4) then
+       print *, "Error@coeff_a_components: expecting width 4 stencil, received width = ", &
+            size(h%stencil)
+       error stop
+    end if
+    if (1 - lbound(h%stencil, 1) /= 2) then
+       print *, "Error@coeff_a_components: expecting centre 2 stencil, received centre = ", &
+            1 - lbound(h%stencil, 1)
+       print *, size(h%stencil), lbound(h%stencil), ubound(h%stencil)
+       error stop
+    end if
+#endif
+
+    select type(deltas => h%stencil)
+    type is(real)
+       hm1 = deltas(-1)
+       h0 = deltas(0)
+       hp1 = deltas(1)
+       hp2 = deltas(2)
+    class default
+       error stop
+    end select
+
+    associate(beta => alpha) ! To match Gamet et al. (1999)
+      numerator = coeff_numerator(hm1, h0, hp2, 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)
+    end associate
+
+  end subroutine coeff_a_components
+  
+  pure real function coeff_numerator(hm1, h0, hp1, hp2, beta)
+    !! Computes the numerator of the coefficient acting on f_{i+1}.
+    !!
+    !! Reduces to (14/3) h^3 when h=const and beta=1/3.
+    
+    real, intent(in) :: hm1 
+    real, intent(in) :: h0 
+    real, intent(in) :: hp1 
+    real, intent(in) :: hp2
+    real, intent(in) :: beta
+    
+    coeff_numerator = h0 * ((hm1 + h0) * (hp1 + hp2) + 2.0 * hp1 * hp2 * beta) ! = (14/3) h^3
+  end function coeff_numerator
+
+  pure real function coeff_numerator_corr(hm1, h0, hp1, hp2, alpha, beta)
+    !! Computes the non-uniform correction to the numerator acting on f_{i+1}.
+    !!
+    !! Reduces to zero when h=const and alpha=beta.
+    
+    real, intent(in) :: hm1 
+    real, intent(in) :: h0 
+    real, intent(in) :: hp1 
+    real, intent(in) :: hp2
+    real, intent(in) :: alpha
+    real, intent(in) :: beta
+
+    coeff_numerator_corr =  hm1 * h0 * hp2 * (beta - alpha) & ! Should cancel for case alpha = beta
+         + (h0**2) * (hp2 - hp1) * beta                     & ! Should cancel for constant h
+         + hm1 * hp1 * (2.0 * hp2 - h0 - hp1) * beta        & ! Should cancel for constant h
+         + (hp1**2) * (3.0 * hp2 - 2.0 * h0 - hp1) * beta   & ! Should cancel for constant h
+         + h0 * (2.0 * hp1 * hp2 * beta - hm1 * (h0 + hp1) * alpha) ! Should cancel for constant h,
+                                                                    ! alpha=beta
+  end function coeff_numerator_corr
+
+  pure real function coeff_denominator(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) :: h0 
+    real, intent(in) :: hp1 
+    real, intent(in) :: hp2 
+
+    coeff_denominator = (3.0 * hp1 * ((h0 + hp1) / 2.0) * hp2)
+  end function coeff_denominator
+
+  pure real function coeff_divisor(hm1, 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))
+  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
new file mode 100644
index 0000000..2e104cc
--- /dev/null
+++ b/src/coeffs/nucfd_coeffs_b.f90
@@ -0,0 +1,172 @@
+! src/coeffs/nucfd_coeff_b.f90
+!
+!! Implements coefficient B of the nucfd_coeffs module.
+!
+! SPDX-License-Identifier: BSD-3-Clause
+
+submodule (nucfd_coeffs) nucfd_coeffs_b
+  !! Submodule defining the coefficient acting on f_{i+1} for compact finite difference schemes on
+  !! non-uniform grids.
+
+  implicit none
+  
+contains
+  
+  module real function coeff_b_points(x)
+    !! Compute the coefficient acting on f_{i-1} of the finite diference given a stencil of points.
+
+    type(nucfd_stencil_points), intent(in) :: x !! Stencil of points for the finite
+                                                !! difference.
+
+    type(nucfd_stencil_deltas) :: h  ! Stencil of grid spacings for the finite difference.
+
+#ifndef NDEBUG
+    if (size(x%stencil) /= 5) then
+       print *, "Error@coeff_b_points: expecting width 5 stencil, received width = ", &
+            size(x%stencil)
+       error stop
+    end if
+    if (1 - lbound(x%stencil, 1) /= 3) then
+       print *, "Error@coeff_b_points: expecting centre 3 stencil, received centre = ", &
+            1 - lbound(x%stencil, 1)
+       print *, size(x%stencil), lbound(x%stencil), ubound(x%stencil)
+       error stop
+    end if
+#endif
+    
+    h = points_to_deltas(x)
+    coeff_b_points = coeff_b_deltas(h)
+    
+  end function coeff_b_points
+
+  module real function coeff_b_deltas(h)
+    !! Compute the coefficient acting on f_{i-1} of the finite diference given a stencil of grid
+    !! spacings.
+
+    type(nucfd_stencil_deltas), intent(in) :: h !! Stencil of grid spacings for the finite difference.
+
+    real :: numerator, numerator_corr, denominator, divisor
+
+#ifndef NDEBUG
+    if (size(h%stencil) /= 4) then
+       print *, "Error@coeff_b_deltas: expecting width 4 stencil, received width = ", &
+            size(h%stencil)
+       error stop
+    end if
+    if (1 - lbound(h%stencil, 1) /= 2) then
+       print *, "Error@coeff_b_deltas: expecting centre 2 stencil, received centre = ", &
+            1 - lbound(h%stencil, 1)
+       print *, size(h%stencil), lbound(h%stencil), ubound(h%stencil)
+       error stop
+    end if
+#endif
+
+    call coeff_b_components(h, numerator, numerator_corr, denominator, divisor)
+    coeff_b_deltas = ((numerator + numerator_corr) / denominator) / divisor
+    
+  end function coeff_b_deltas
+
+  module subroutine coeff_b_components(h, numerator, numerator_corr, denominator, divisor)
+
+    type(nucfd_stencil_deltas), intent(in) :: h !! Stencil of grid spacings for the finite difference.
+    real, intent(out) :: numerator
+    real, intent(out) :: numerator_corr
+    real, intent(out) :: denominator
+    real, intent(out) :: divisor
+
+    real :: hm1, h0, hp1, hp2 ! Grid deltas at i -2, -1, 0, +1, +2
+
+#ifndef NDEBUG
+    if (size(h%stencil) /= 4) then
+       print *, "Error@coeff_b_components: expecting width 4 stencil, received width = ", &
+            size(h%stencil)
+       error stop
+    end if
+    if (1 - lbound(h%stencil, 1) /= 2) then
+       print *, "Error@coeff_b_components: expecting centre 2 stencil, received centre = ", &
+            1 - lbound(h%stencil, 1)
+       print *, size(h%stencil), lbound(h%stencil), ubound(h%stencil)
+       error stop
+    end if
+#endif
+
+    select type(deltas => h%stencil)
+    type is(real)
+       hm1 = deltas(-1)
+       h0 = deltas(0)
+       hp1 = deltas(1)
+       hp2 = deltas(2)
+    class default
+       error stop
+    end select
+
+    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)
+    end associate
+
+  end subroutine coeff_b_components
+
+  pure real function coeff_numerator(hm1, h0, hp1, hp2, alpha)
+    !! Computes the numerator of the coefficient acting on f_{i-1}.
+    !!
+    !! Reduces to -(14/3) h^3 when h=const and alpha=beta.
+    
+    real, intent(in) :: hm1 
+    real, intent(in) :: h0 
+    real, intent(in) :: hp1 
+    real, intent(in) :: hp2
+    real, intent(in) :: alpha
+    
+    coeff_numerator = -hp1 * ((hm1 + h0) * (hp1 + hp2) + 2.0 * hm1 * h0 * alpha) ! = -(14/3) h^3
+  end function coeff_numerator
+
+  pure real function coeff_numerator_corr(hm1, h0, hp1, hp2, alpha, beta)
+    !! Computes the non-uniform correction to the numerator acting on f_{i-1}.
+    !!
+    !! Reduces to zero when h=const and alpha=beta.
+    
+    real, intent(in) :: hm1 
+    real, intent(in) :: h0 
+    real, intent(in) :: hp1 
+    real, intent(in) :: hp2
+    real, intent(in) :: alpha
+    real, intent(in) :: beta
+
+    coeff_numerator_corr = -hm1 * hp1 * hp2 * (alpha - beta) & ! Should cancel for case alpha = beta
+         - (hp1**2) * (hm1 - h0) * alpha &    ! Should cancel for case h = const
+         - h0 * hp2 * (2.0 * hm1 - h0 - hp1) * alpha &      ! Should cancel for case h = const
+         - (h0**2) * (3.0 * hm1 - h0 - 2.0 * hp1) * alpha & ! Should cancel for case h = const
+         - hp1 * (2.0 * hm1 * h0 * alpha - (h0 + hp1) * hp2 * beta) ! Should cancel for constant h &
+                                                                    ! alpha = beta
+  end function coeff_numerator_corr
+
+  pure real function coeff_denominator(hm1, h0, hp1)
+    !! 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
+    
+    coeff_denominator = 3.0 * hm1 * h0 * ((h0 + hp1) / 2.0)
+  end function coeff_denominator
+
+  pure real function coeff_divisor(h0, hp1, hp2)
+    !! 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) :: h0
+    real, intent(in) :: hp1
+    real, intent(in) :: hp2
+
+    coeff_divisor = 2.0 * (h0 + hp1 + hp2) / 3.0
+  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
new file mode 100644
index 0000000..256a716
--- /dev/null
+++ b/src/coeffs/nucfd_coeffs_c.f90
@@ -0,0 +1,161 @@
+! src/coeffs/nucfd_coeff_c.f90
+!
+!! Implements coefficient C of the nucfd_coeffs module.
+!
+! SPDX-License-Identifier: BSD-3-Clause
+
+submodule (nucfd_coeffs) nucfd_coeffs_c
+  !! Submodule defining the coefficient acting on f_{i+2} for compact finite difference schemes on
+  !! non-uniform grids.
+
+  implicit none
+
+contains
+
+  module real function coeff_c_points(x)
+    !! Compute the coefficient acting on f_{i+2} of the finite diference given a stencil of points.
+
+    type(nucfd_stencil_points), intent(in) :: x !! Stencil of points for the finite
+                                                !! difference.
+
+    type(nucfd_stencil_deltas) :: h  ! Stencil of grid spacings for the finite difference.
+
+#ifndef NDEBUG
+    if (size(x%stencil) /= 5) then
+       print *, "Error@coeff_c_points: expecting width 5 stencil, received width = ", &
+            size(x%stencil)
+       error stop
+    end if
+    if (1 - lbound(x%stencil, 1) /= 3) then
+       print *, "Error@coeff_c_points: expecting centre 3 stencil, received centre = ", &
+            1 - lbound(x%stencil, 1)
+       print *, size(x%stencil), lbound(x%stencil), ubound(x%stencil)
+       error stop
+    end if
+#endif
+    
+    h = points_to_deltas(x)
+    coeff_c_points = coeff_c_deltas(h)
+    
+  end function coeff_c_points
+
+  module real function coeff_c_deltas(h)
+    !! Compute the coefficient acting on f_{i+2} of the finite diference given a stencil of grid
+    !! spacings.
+
+    type(nucfd_stencil_deltas), intent(in) :: h !! Stencil of grid spacings for the finite difference.
+
+    real :: numerator, denominator, divisor
+
+#ifndef NDEBUG
+    if (size(h%stencil) /= 4) then
+       print *, "Error@coeff_c_deltas: expecting width 4 stencil, received width = ", &
+            size(h%stencil)
+       error stop
+    end if
+    if (1 - lbound(h%stencil, 1) /= 2) then
+       print *, "Error@coeff_c_deltas: expecting centre 2 stencil, received centre = ", &
+            1 - lbound(h%stencil, 1)
+       print *, size(h%stencil), lbound(h%stencil), ubound(h%stencil)
+       error stop
+    end if
+#endif
+    
+    call coeff_c_components(h, numerator, denominator, divisor)
+    coeff_c_deltas = (numerator / denominator) / divisor
+    
+  end function coeff_c_deltas
+
+  module subroutine coeff_c_components(h, numerator, denominator, divisor)
+    !! Compute the components of the coefficient acting on f_{i+2} of the finite diference given a
+    !! stencil of grid spacings.
+    !!
+    !! For uniform grids the numerator should reduce to (2/3) h^3, the denominator to 6h^3 (for a
+    !! coefficient value of 1/9) and the finite difference divisor to 4h.
+
+    type(nucfd_stencil_deltas), intent(in) :: h !! Stencil of grid spacings for the finite difference.
+    real, intent(out) :: numerator   !! The numerator of the coefficient
+    real, intent(out) :: denominator !! The denominator of the coefficient
+    real, intent(out) :: divisor     !! The finite-difference divisor of the coefficient
+    
+    real :: hm1, h0, hp1, hp2 ! Grid deltas at i -1, 0, +1, +2
+
+#ifndef NDEBUG
+    if (size(h%stencil) /= 4) then
+       print *, "Error@coeff_c_components: expecting width 4 stencil, received width = ", &
+            size(h%stencil)
+       error stop
+    end if
+    if (1 - lbound(h%stencil, 1) /= 2) then
+       print *, "Error@coeff_c_components: expecting centre 2 stencil, received centre = ", &
+            1 - lbound(h%stencil, 1)
+       print *, size(h%stencil), lbound(h%stencil), ubound(h%stencil)
+       error stop
+    end if
+#endif
+
+    select type(deltas => h%stencil)
+    type is(real)
+       hm1 = deltas(-1)
+       h0 = deltas(0)
+       hp1 = deltas(1)
+       hp2 = deltas(2)
+    class default
+       error stop
+    end select
+
+    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)
+    end associate
+    
+  end subroutine coeff_c_components
+    
+  pure real function coeff_numerator(hm1, h0, hp1, alpha, beta)
+    !! Computes the numerator of the coefficient acting on f_{i+2}.
+    !!
+    !! Reduces to (2/3) h^3 when h=const and alpha=beta=1/3.
+
+    real, intent(in) :: hm1
+    real, intent(in) :: h0
+    real, intent(in) :: hp1
+    real, intent(in) :: alpha
+    real, intent(in) :: beta
+    
+    coeff_numerator = -h0 * hp1 * (hm1 + h0) &
+         + hm1 * h0 * (h0 + hp1) * alpha &
+         + hp1 * (h0 * (hm1 + h0) + hp1 * (hm1 + 2.0 * h0 + hp1)) * beta
+    
+  end function coeff_numerator
+  
+  pure real function coeff_denominator(hm1, 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)
+    
+  end function coeff_denominator
+
+  pure real function coeff_divisor(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) :: h0
+    real, intent(in) :: hp1
+    real, intent(in) :: hp2
+
+    coeff_divisor = 4.0 * ((h0 + hp1 + hp2) / 3.0) 
+  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
new file mode 100644
index 0000000..9dd0f80
--- /dev/null
+++ b/src/coeffs/nucfd_coeffs_d.f90
@@ -0,0 +1,161 @@
+! src/coeffs/nucfd_coeff_d.f90
+!
+!! Implements coefficient D of the nucfd_coeffs module.
+!
+! SPDX-License-Identifier: BSD-3-Clause
+
+submodule (nucfd_coeffs) nucfd_coeffs_d
+  !! Submodule defining the coefficient acting on f_{i+2} for compact finite difference schemes on
+  !! non-uniform grids.
+
+  implicit none
+
+contains
+  
+  module real function coeff_d_points(x)
+    !! Compute the coefficient acting on f_{i-2} of the finite diference given a stencil of points.
+
+    type(nucfd_stencil_points), intent(in) :: x !! Stencil of points for the finite
+                                                !! difference.
+
+    type(nucfd_stencil_deltas) :: h  ! Stencil of grid spacings for the finite difference.
+
+#ifndef NDEBUG
+    if (size(x%stencil) /= 5) then
+       print *, "Error@coeff_d_points: expecting width 5 stencil, received width = ", &
+            size(x%stencil)
+       error stop
+    end if
+    if (1 - lbound(x%stencil, 1) /= 3) then
+       print *, "Error@coeff_d_points: expecting centre 3 stencil, received centre = ", &
+            1 - lbound(x%stencil, 1)
+       print *, size(x%stencil), lbound(x%stencil), ubound(x%stencil)
+       error stop
+    end if
+#endif
+    
+    h = points_to_deltas(x)
+    coeff_d_points = coeff_d_deltas(h)
+    
+  end function coeff_d_points
+
+  module real function coeff_d_deltas(h)
+    !! Compute the coefficient acting on f_{i-2} of the finite diference given a stencil of grid
+    !! spacings.
+
+    type(nucfd_stencil_deltas), intent(in) :: h !! Stencil of grid spacings for the finite difference.
+
+    real :: numerator, denominator, divisor
+
+#ifndef NDEBUG
+    if (size(h%stencil) /= 4) then
+       print *, "Error@coeff_d_deltas: expecting width 4 stencil, received width = ", &
+            size(h%stencil)
+       error stop
+    end if
+    if (1 - lbound(h%stencil, 1) /= 2) then
+       print *, "Error@coeff_d_deltas: expecting centre 2 stencil, received centre = ", &
+            1 - lbound(h%stencil, 1)
+       print *, size(h%stencil), lbound(h%stencil), ubound(h%stencil)
+       error stop
+    end if
+#endif
+    
+    call coeff_d_components(h, numerator, denominator, divisor)
+    coeff_d_deltas = (numerator / denominator) / divisor
+    
+  end function coeff_d_deltas
+
+  module subroutine coeff_d_components(h, numerator, denominator, divisor)
+    !! Compute the components of the coefficient acting on f_{i-2} of the finite diference given a
+    !! stencil of grid spacings.
+    !!
+    !! For uniform grids the numerator should reduce to -(2/3) h^3, the denominator to 6h^3 (for a
+    !! coefficient value of -1/9) and the finite difference divisor to 4h.
+
+    type(nucfd_stencil_deltas), intent(in) :: h !! Stencil of grid spacings for the finite difference.
+    real, intent(out) :: numerator   !! The numerator of the coefficient
+    real, intent(out) :: denominator !! The denominator of the coefficient
+    real, intent(out) :: divisor     !! The finite-difference divisor of the coefficient
+    
+    real :: hm1, h0, hp1, hp2 ! Grid deltas at i -1, 0, +1, +2
+
+#ifndef NDEBUG
+    if (size(h%stencil) /= 4) then
+       print *, "Error@coeff_d_components: expecting width 4 stencil, received width = ", &
+            size(h%stencil)
+       error stop
+    end if
+    if (1 - lbound(h%stencil, 1) /= 2) then
+       print *, "Error@coeff_d_components: expecting centre 2 stencil, received centre = ", &
+            1 - lbound(h%stencil, 1)
+       print *, size(h%stencil), lbound(h%stencil), ubound(h%stencil)
+       error stop
+    end if
+#endif
+
+    select type(deltas => h%stencil)
+    type is(real)
+       hm1 = deltas(-1)
+       h0 = deltas(0)
+       hp1 = deltas(1)
+       hp2 = deltas(2)
+    class default
+       error stop
+    end select
+
+    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)
+    end associate
+    
+  end subroutine coeff_d_components
+
+  pure real function coeff_numerator(h0, hp1, hp2, alpha, beta)
+    !! Computes the numerator of the coefficient acting on f_{i-2}.
+    !!
+    !! Reduces to -(2/3) h^3 when h=const and alpha=beta=1/3.
+
+    real, intent(in) :: h0
+    real, intent(in) :: hp1
+    real, intent(in) :: hp2
+    real, intent(in) :: alpha
+    real, intent(in) :: beta
+    
+    coeff_numerator = h0 * hp1 * (hp1 + hp2) &
+         - h0 * (h0 * (h0 + 2.0 * hp1 + hp2) + hp1 * (hp1 + hp2)) * alpha &
+         - hp1 * hp2 * (h0 + hp1) * beta
+    
+  end function coeff_numerator
+  
+  pure real function coeff_denominator(hm1, h0, hp1, hp2)
+    !! 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
+    !! -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 * hm1 * (hp1 + hp2) &
+         * ((hm1 + h0 + hp1 + hp2) / 4.0)
+    
+  end function coeff_denominator
+
+  pure real function coeff_divisor(hm1, h0, hp1)
+    !! 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
+
+    coeff_divisor = 4.0 * ((hm1 + h0 + hp1) / 3.0) 
+  end function coeff_divisor
+
+end submodule nucfd_coeffs_d
diff --git a/src/coeffs/nucfd_coeffs_e.f90 b/src/coeffs/nucfd_coeffs_e.f90
new file mode 100644
index 0000000..824ffd9
--- /dev/null
+++ b/src/coeffs/nucfd_coeffs_e.f90
@@ -0,0 +1,46 @@
+! src/coeffs/nucfd_coeff_e.f90
+!
+!! Implements coefficient E of the nucfd_coeffs module.
+!
+! SPDX-License-Identifier: BSD-3-Clause
+
+submodule (nucfd_coeffs) nucfd_coeffs_e
+  !! Submodule defining the coefficient acting on f_i for compact finite difference schemes on
+  !! non-uniform grids.
+
+  implicit none
+
+contains
+
+  module real function coeff_e(x)
+    !! Compute the coefficient acting on f_{i} of the finite diference given a stencil of points.
+
+    type(nucfd_stencil_points), intent(in) :: x !! Stencil of points for the finite
+                                                !! difference.
+
+    real :: a, b, c, d ! The neighbouring coefficients
+
+#ifndef NDEBUG
+    if (size(x%stencil) /= 5) then
+       print *, "Error@coeff_e: expecting width 5 stencil, received width = ", &
+            size(x%stencil)
+       error stop
+    end if
+    if (1 - lbound(x%stencil, 1) /= 3) then
+       print *, "Error@coeff_e: expecting centre 3 stencil, received centre = ", &
+            1 - lbound(x%stencil, 1)
+       print *, size(x%stencil), lbound(x%stencil), ubound(x%stencil)
+       error stop
+    end if
+#endif
+
+    a = coeff_a(x)
+    b = coeff_b(x)
+    c = coeff_c(x)
+    d = coeff_d(x)
+
+    coeff_e = -(a + b + c + d)
+    
+  end function coeff_e
+  
+end submodule nucfd_coeffs_e
diff --git a/src/coeffs/nucfd_coeffs_mod.f90 b/src/coeffs/nucfd_coeffs_mod.f90
new file mode 100644
index 0000000..6c57d55
--- /dev/null
+++ b/src/coeffs/nucfd_coeffs_mod.f90
@@ -0,0 +1,127 @@
+! src/coeffs/nucfd_coeff_mod.f90
+!
+!! Defines the nucfd_coeffs module.
+!
+! SPDX-License-Identifier: BSD-3-Clause
+
+module nucfd_coeffs
+  !! Module defining the coefficients for compact finite difference schemes on non-uniform grids.
+
+  use nucfd_types
+  
+  implicit none
+
+  private
+  public :: coeff_a, coeff_a_components
+  public :: coeff_b, coeff_b_components
+  public :: coeff_c, coeff_c_components
+  public :: coeff_d, coeff_d_components
+  public :: coeff_e
+
+  interface coeff_a
+     !! Compute the coefficient acting on f_{i+1} of the finite diference stencil.
+     module procedure coeff_a_points
+     module procedure coeff_a_deltas
+  end interface coeff_a
+  
+  interface coeff_b
+     !! Compute the coefficient acting on f_{i-1} of the finite diference stencil.
+     module procedure coeff_b_points
+     module procedure coeff_b_deltas
+  end interface coeff_b
+
+  interface coeff_c
+     !! Compute the coefficient acting on f_{i+2} of the finite diference stencil.
+     module procedure coeff_c_points
+     module procedure coeff_c_deltas
+  end interface coeff_c
+
+  interface coeff_d
+     !! Compute the coefficient acting on f_{i-2} of the finite diference stencil.
+     module procedure coeff_d_points
+     module procedure coeff_d_deltas
+  end interface coeff_d
+
+  interface
+     module real function coeff_a_points(x)
+       type(nucfd_stencil_points), intent(in) :: x !! Stencil of points for the finite
+                                                   !! difference.
+     end function
+     module real function coeff_a_deltas(h)
+       type(nucfd_stencil_deltas), intent(in) :: h !! Stencil of grid spacings for the finite
+                                                   !! difference.
+     end function
+     module subroutine coeff_a_components(h, numerator, numerator_corr, denominator, divisor)
+       type(nucfd_stencil_deltas), intent(in) :: h !! Stencil of grid spacings for the finite difference.
+       real, intent(out) :: numerator
+       real, intent(out) :: numerator_corr
+       real, intent(out) :: denominator
+       real, intent(out) :: divisor
+     end subroutine coeff_a_components
+
+     module real function coeff_b_points(x)
+       type(nucfd_stencil_points), intent(in) :: x !! Stencil of points for the finite
+                                                   !! difference.
+     end function
+     module real function coeff_b_deltas(h)
+       type(nucfd_stencil_deltas), intent(in) :: h !! Stencil of grid spacings for the finite
+                                                   !! difference.
+     end function
+     module subroutine coeff_b_components(h, numerator, numerator_corr, denominator, divisor)
+       type(nucfd_stencil_deltas), intent(in) :: h !! Stencil of grid spacings for the finite difference.
+       real, intent(out) :: numerator
+       real, intent(out) :: numerator_corr
+       real, intent(out) :: denominator
+       real, intent(out) :: divisor
+     end subroutine coeff_b_components
+
+     module real function coeff_c_points(x)
+       type(nucfd_stencil_points), intent(in) :: x !! Stencil of points for the finite
+                                                   !! difference.
+     end function
+     module real function coeff_c_deltas(h)
+       type(nucfd_stencil_deltas), intent(in) :: h !! Stencil of grid spacings for the finite
+                                                   !! difference.
+     end function
+     module subroutine coeff_c_components(h, numerator, denominator, divisor)
+       !! Compute the components of the coefficient acting on f_{i+2} of the finite diference given a
+       !! stencil of grid spacings.
+       !!
+       !! For uniform grids the numerator should reduce to (2/3) h^3, the denominator to 6h^3 (for a
+       !! coefficient value of (1/9) h^3) and the finite difference divisor to 4h.
+       type(nucfd_stencil_deltas), intent(in) :: h !! Stencil of grid spacings for the finite difference.
+       real, intent(out) :: numerator   !! The numerator of the coefficient
+       real, intent(out) :: denominator !! The denominator of the coefficient
+       real, intent(out) :: divisor     !! The finite-difference divisor of the coefficient
+     end subroutine coeff_c_components
+     
+     module real function coeff_d_points(x)
+       type(nucfd_stencil_points), intent(in) :: x !! Stencil of points for the finite
+                                                   !! difference.
+     end function
+     module real function coeff_d_deltas(h)
+       type(nucfd_stencil_deltas), intent(in) :: h !! Stencil of grid spacings for the finite
+                                                   !! difference.
+     end function
+     module subroutine coeff_d_components(h, numerator, denominator, divisor)
+       !! Compute the components of the coefficient acting on f_{i-2} of the finite diference given a
+       !! stencil of grid spacings.
+       !!
+       !! For uniform grids the numerator should reduce to -(2/3) h^3, the denominator to 6h^3 (for a
+       !! coefficient value of -1/9) and the finite difference divisor to 4h.
+       type(nucfd_stencil_deltas), intent(in) :: h !! Stencil of grid spacings for the finite difference.
+       real, intent(out) :: numerator   !! The numerator of the coefficient
+       real, intent(out) :: denominator !! The denominator of the coefficient
+       real, intent(out) :: divisor     !! The finite-difference divisor of the coefficient
+     end subroutine coeff_d_components
+  
+     module real function coeff_e(x)
+       type(nucfd_stencil_points), intent(in) :: x !! Stencil of points for the finite
+                                                   !! difference.
+     end function
+  end interface
+  
+  real, parameter, public :: alpha = 1.0 / 3.0 !! Off-diagonal coefficient for first derivative
+                                               !! system.
+
+end module nucfd_coeffs
diff --git a/src/nucfd_coeffs_mod.f90 b/src/nucfd_coeffs_mod.f90
deleted file mode 100644
index 82135f5..0000000
--- a/src/nucfd_coeffs_mod.f90
+++ /dev/null
@@ -1,13 +0,0 @@
-module nucfd_coeffs
-  !! Module defining the coefficients for compact finite difference schemes on non-uniform grids.
-  !!
-  !! SPDX-License-Identifier: BSD-3-Clause
-
-  implicit none
-
-  private
-
-  real, parameter, public :: alpha = 1.0 / 3.0 !! Off-diagonal coefficient for first derivative
-                                               !! system.
-
-end module nucfd_coeffs
diff --git a/src/nucfd_deriv_mod.f90 b/src/nucfd_deriv_mod.f90
new file mode 100644
index 0000000..f9a946b
--- /dev/null
+++ b/src/nucfd_deriv_mod.f90
@@ -0,0 +1,64 @@
+! src/nucfd_deriv_mod.f90
+!
+!! Defines the nucfd_deriv module.
+!
+! SPDX-License-Identifier: BSD-3-Clause
+
+module nucfd_deriv
+  !! Module defining the coefficients for compact finite difference schemes on non-uniform grids.
+
+  use nucfd_types
+  use nucfd_coeffs
+  
+  implicit none
+
+  private
+  public :: deriv_rhs
+
+contains
+  
+  subroutine deriv_rhs(f, stencil, x, dfdx)
+    !! Compute the RHS of the derivative of a function f based on a stencil of indices.
+
+    real, dimension(:), intent(in) :: f
+    type(nucfd_index_stencil), intent(in) :: stencil
+    real, dimension(:), intent(in) :: x
+    real, intent(out) :: dfdx
+
+    type(nucfd_stencil_points) :: stencil_coordinates
+
+    real :: a, b, c, d, e
+    
+    call create_stencil(5, 3, stencil_coordinates)
+
+    select type(indices => stencil%stencil)
+    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))
+       class default
+          print *, "Error: Coordinate stencil is misallocated!"
+          error stop
+       end select
+
+       a = coeff_a(stencil_coordinates)
+       b = coeff_b(stencil_coordinates)
+       c = coeff_c(stencil_coordinates)
+       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))
+    class default
+       print *, "Error: Index stencil is misallocated!"
+       error stop
+    end select
+    
+  end subroutine deriv_rhs
+  
+end module nucfd_deriv
diff --git a/src/nucfd_trid_solver_mod.f90 b/src/nucfd_trid_solver_mod.f90
index 892f067..97e4c40 100644
--- a/src/nucfd_trid_solver_mod.f90
+++ b/src/nucfd_trid_solver_mod.f90
@@ -1,9 +1,13 @@
+! src/nucfd_trid_solver_mod.f90
+!
+!! Defines the nucfd_trid_solver module.
+!
+! SPDX-License-Identifier: BSD-3-Clause
+
 module nucfd_trid_solver
   !! Module implementing a simple tridiagonal solver based on the algorithm given at
   !! https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm including support for cyclic
   !! systems.
-  !!
-  !! SPDX-License-Identifier: BSD-3-Clause
 
   implicit none
 
diff --git a/src/nucfd_types_mod.f90 b/src/nucfd_types_mod.f90
new file mode 100644
index 0000000..1c2bace
--- /dev/null
+++ b/src/nucfd_types_mod.f90
@@ -0,0 +1,112 @@
+! src/nucfd_types_mod.f90
+!
+!! Defines the nucfd_types module.
+!
+! SPDX-License-Identifier: BSD-3-Clause
+
+module nucfd_types
+  !! Module defining the custom types used by NuCFD.
+
+  implicit none
+
+  type nucfd_stencil
+     !! Type representing a stencil.
+     class(*), dimension(:), allocatable :: stencil !! Stencil data
+   contains
+     final :: free_stencil
+  end type nucfd_stencil
+  
+  type, extends(nucfd_stencil) :: nucfd_stencil_points
+     !! Type representing the point coordinates of a stencil.
+  end type nucfd_stencil_points
+
+  type, extends(nucfd_stencil) :: nucfd_stencil_deltas
+     !! Type representing the grid spacings of a stencil. Following Gamet et al. (1999) these are
+     !! defined as h_i = x_i - x_{i-1}.
+  end type nucfd_stencil_deltas
+
+  type, extends(nucfd_stencil) :: nucfd_index_stencil
+     !! Type representing a stencil's indices.
+  end type nucfd_index_stencil
+
+contains
+
+  subroutine create_stencil(width, centre, stencil)
+    !! Subroutine to create a stencil, specifying the width and centre. The underlying array is
+    !! allocated with bounds (1-centre):(1-centre)+(width-1) so that accesses are by offset from the
+    !! central position.
+
+    integer, intent(in) :: width                 !! Specifies the width of the stencil
+    integer, intent(in) :: centre                !! Specifies the index the stencil should be centred about
+    class(nucfd_stencil), intent(out) :: stencil !! The stencil object.
+
+    integer :: lb, ub
+
+    lb = 1 - centre
+    ub = (1 - centre) + (width - 1)
+    select type(stencil)
+    type is(nucfd_stencil_points)
+       allocate(real::stencil%stencil(lb:ub))
+    type is(nucfd_stencil_deltas)
+       allocate(real::stencil%stencil(lb:ub))
+    type is(nucfd_index_stencil)
+       allocate(integer::stencil%stencil(lb:ub))
+    class default
+       print *, "Error: trying to create unknown stencil type"
+       error stop
+    end select
+    
+  end subroutine create_stencil
+  
+  function points_to_deltas(x) result(h)
+    !! Helper function converting a stencil's point representation to a grid spacing representation.
+
+    type(nucfd_stencil_points), intent(in) :: x !! The stencil's point representation.
+    type(nucfd_stencil_deltas) :: h             !! The stencil's grid spacing representation.
+
+    integer :: i
+    real :: xm1, x0 ! Start and end points of grid spacing at i.
+
+    integer :: ub, lb, width, centre ! Upper bound, lower bound, width and centre of grid stencil.
+
+    lb = lbound(x%stencil, 1) + 1
+    ub = ubound(x%stencil, 1)
+    width = (ub - lb) + 1
+    centre = 1 - lb
+    call create_stencil(width, centre, h)
+
+    associate(points => x%stencil, &
+         deltas => h%stencil)
+      select type(points)
+      type is(real)
+         select type(deltas)
+         type is(real)
+            do i = lb, ub
+               xm1 = points(i - 1)
+               x0 = points(i)
+
+               deltas(i) = x0 - xm1
+            end do
+         class default
+            print *, "Error: Difference stencil is misallocated!"
+         end select
+      class default
+         print *, "Error: Coordinate stencil is misallocated!"
+         error stop
+      end select
+    end associate
+    
+   end function points_to_deltas
+
+   subroutine free_stencil(stencil)
+     !! Frees a stencil object
+     
+     type(nucfd_stencil) :: stencil !! The stencil object to be freed.
+
+     if (allocated(stencil%stencil)) then
+        deallocate(stencil%stencil)
+     end if
+     
+   end subroutine free_stencil
+   
+end module nucfd_types
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
new file mode 100644
index 0000000..d2342c1
--- /dev/null
+++ b/tests/CMakeLists.txt
@@ -0,0 +1,37 @@
+## tests/CMakeLists.txt
+#
+## Description
+#
+# CMakeLists file controlling the tests for the NuCFD project.
+#
+## LICENSE
+#
+# SPDX-License-Identifier: BSD-3-Clause
+#
+
+add_library(nucfd_test_flags INTERFACE)
+target_link_libraries(nucfd_test_flags INTERFACE
+  nucfd_compiler_flags
+  nucfd_debug_flags
+  nucfd_strict_flags)
+
+add_library(nucfd_test_framework 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} nucfd_test_flags)
+  add_test(NAME ${suite}:${test_name} COMMAND ${test_name})
+endfunction()
+
+function(add_test_libs suite test_name lib)
+  target_link_libraries(${test_name} ${lib})
+endfunction()
+
+# Tridiagonal solver test suite
+add_subdirectory(tridsolver)
+
+# Finite difference test suite
+add_subdirectory(fd-schemes)
diff --git a/tests/fd-schemes/CMakeLists.txt b/tests/fd-schemes/CMakeLists.txt
new file mode 100644
index 0000000..27bea18
--- /dev/null
+++ b/tests/fd-schemes/CMakeLists.txt
@@ -0,0 +1,13 @@
+## tests/fd-schemes/CMakeLists.txt
+#
+## Description
+#
+# CMakeLists file for the finite difference schemes test suite.
+#
+## LICENSE
+#
+# SPDX-License-Identifier: BSD-3-Clause
+#
+
+add_subdirectory(verify_coeffs)
+add_subdirectory(differentiation_rules)
diff --git a/tests/fd-schemes/differentiation_rules/CMakeLists.txt b/tests/fd-schemes/differentiation_rules/CMakeLists.txt
new file mode 100644
index 0000000..3ffc88a
--- /dev/null
+++ b/tests/fd-schemes/differentiation_rules/CMakeLists.txt
@@ -0,0 +1,14 @@
+## tests/fd-schemes/differentiation_rules/CMakeLists.txt
+#
+## Description
+#
+# CMakeLists file for the finite difference schemes differentiation rules test sub-suite.
+#
+## LICENSE
+#
+# SPDX-License-Identifier: BSD-3-Clause
+#
+
+define_test(fd-schemes constant_function)
+define_test(fd-schemes linear_rising_function)
+define_test(fd-schemes linear_falling_function)
diff --git a/tests/fd-schemes/differentiation_rules/constant_function.f90 b/tests/fd-schemes/differentiation_rules/constant_function.f90
new file mode 100644
index 0000000..06ea223
--- /dev/null
+++ b/tests/fd-schemes/differentiation_rules/constant_function.f90
@@ -0,0 +1,77 @@
+! tests/fd-schemes/differentiation_rules/constant_function.f90
+!
+!! Part of the fd-schemes test suite.
+!
+! SPDX-License-Identifier: BSD-3-Clause
+
+program constant_function
+  !! Tests that rules of differentiation are respected for constant functions.
+
+  use nucfd_types
+  use nucfd_deriv
+  
+  use nucfd_tests
+  
+  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
+  
+  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
+
+  f(:) = 1.0
+  call deriv_rhs(f, stencil, x, dfdx)
+  call deriv_rhs(f + 1.0, stencil, x, dgdx)
+  call test_report("Shifted derivative +, const f", check_scalar(dgdx, dfdx))
+  call deriv_rhs(f - 1.0, stencil, x, dgdx)
+  call test_report("Shifted derivative -, const f", check_scalar(dgdx, dfdx))
+  call deriv_rhs(f, stencil, x, dgdx)
+  call test_report("Shifted derivative 0, const f", check_scalar(dgdx, dfdx))
+  call deriv_rhs(2.0 * f, stencil, x, dgdx)
+  call test_report("Scaled derivative 2x, const f", check_scalar(dgdx, 2.0 * dfdx))
+  call deriv_rhs(-f, stencil, x, dgdx)
+  call test_report("Scaled derivative -1, const f", check_scalar(dgdx, -dfdx))
+  
+  deallocate(x)
+  deallocate(f)
+  
+  call finalise_suite()
+
+end program constant_function
diff --git a/tests/fd-schemes/differentiation_rules/linear_falling_function.f90 b/tests/fd-schemes/differentiation_rules/linear_falling_function.f90
new file mode 100644
index 0000000..eab0827
--- /dev/null
+++ b/tests/fd-schemes/differentiation_rules/linear_falling_function.f90
@@ -0,0 +1,82 @@
+! tests/fd-schemes/differentiation_rules/linear_falling_function.f90
+!
+!! Part of the fd-schemes test suite.
+!
+! SPDX-License-Identifier: BSD-3-Clause
+
+program linear_falling_function
+  !! Tests that rules of differentiation are respected for linearly decreasing functions.
+
+  use nucfd_types
+  use nucfd_deriv
+  
+  use nucfd_tests
+  
+  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
+  
+  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
+
+  f(1) = 0.0
+  dfdx = 1.0
+  do i = 2, n
+     f(i) = f(i - 1) - h * dfdx
+  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 linear_falling_function
diff --git a/tests/fd-schemes/differentiation_rules/linear_rising_function.f90 b/tests/fd-schemes/differentiation_rules/linear_rising_function.f90
new file mode 100644
index 0000000..4400537
--- /dev/null
+++ b/tests/fd-schemes/differentiation_rules/linear_rising_function.f90
@@ -0,0 +1,81 @@
+! tests/fd-schemes/differentiation_rules/linear_rising_function.f90
+!
+!! Part of the fd-schemes test suite.
+!
+! SPDX-License-Identifier: BSD-3-Clause
+
+program linear_rising_function
+  !! Tests that rules of differentiation are respected for linearly increasing functions.
+
+  use nucfd_types
+  use nucfd_deriv
+  
+  use nucfd_tests
+  
+  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
+  
+  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
+
+  f(1) = 0.0
+  dfdx = 1.0
+  do i = 2, n
+     f(i) = f(i - 1) + h * dfdx
+  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 linear_rising_function
diff --git a/tests/fd-schemes/verify_coeffs/CMakeLists.txt b/tests/fd-schemes/verify_coeffs/CMakeLists.txt
new file mode 100644
index 0000000..f95f68c
--- /dev/null
+++ b/tests/fd-schemes/verify_coeffs/CMakeLists.txt
@@ -0,0 +1,16 @@
+## tests/fd-schemes/verify_coeffs/CMakeLists.txt
+#
+## Description
+#
+# CMakeLists file for the finite difference schemes differentiation rules test sub-suite.
+#
+## LICENSE
+#
+# SPDX-License-Identifier: BSD-3-Clause
+#
+
+define_test(fd-schemes verify_coeff_a)
+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)
diff --git a/tests/fd-schemes/verify_coeffs/verify_coeff_a.f90 b/tests/fd-schemes/verify_coeffs/verify_coeff_a.f90
new file mode 100644
index 0000000..9c3e5b4
--- /dev/null
+++ b/tests/fd-schemes/verify_coeffs/verify_coeff_a.f90
@@ -0,0 +1,61 @@
+! tests/fd-schemes/verify_coeffs/verify_coeff_a.f90
+!
+!! Part of the fd-schemes test suite.
+!
+! SPDX-License-Identifier: BSD-3-Clause
+
+program verify_coeff_b
+  !! Tests the computation of finite difference coefficient for non-uniform grids acting on f_{i+1}.
+
+  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
+  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
+  
+  call initialise_suite("Verify coefficient A")
+
+  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
+
+  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)))
+
+  call finalise_suite()
+
+end program verify_coeff_b
diff --git a/tests/fd-schemes/verify_coeffs/verify_coeff_b.f90 b/tests/fd-schemes/verify_coeffs/verify_coeff_b.f90
new file mode 100644
index 0000000..dcd013f
--- /dev/null
+++ b/tests/fd-schemes/verify_coeffs/verify_coeff_b.f90
@@ -0,0 +1,61 @@
+! tests/fd-schemes/verify_coeffs/verify_coeff_b.f90
+!
+!! Part of the fd-schemes test suite.
+!
+! SPDX-License-Identifier: BSD-3-Clause
+
+program verify_coeff_b
+  !! Tests the computation of finite difference coefficient for non-uniform grids acting on f_{i-1}.
+
+  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 :: b
+  real, parameter :: bref = -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
+  
+  call initialise_suite("Verify coefficient B")
+
+  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
+
+  call coeff_b_components(points_to_deltas(stencil), numerator, numerator_corr, denominator, divisor)
+  call test_report("Coefficient B numerator", check_scalar(numerator, numerator_f1ref * (h**3)))
+  call test_report("Coefficient B numerator correction", check_scalar(numerator_corr, numerator_corr_f1ref * (h**3)))
+  call test_report("Coefficient B denominator", check_scalar(denominator, denominator_f1ref * (h**3)))
+  call test_report("Coefficient B divisor", check_scalar(divisor, divisor_f1ref * h))
+  b = coeff_b(stencil)
+  call test_report("Coefficient B", check_scalar(b, bref / (2.0 * h)))
+
+  call finalise_suite()
+
+end program verify_coeff_b
diff --git a/tests/fd-schemes/verify_coeffs/verify_coeff_c.f90 b/tests/fd-schemes/verify_coeffs/verify_coeff_c.f90
new file mode 100644
index 0000000..16a8787
--- /dev/null
+++ b/tests/fd-schemes/verify_coeffs/verify_coeff_c.f90
@@ -0,0 +1,59 @@
+! tests/fd-schemes/verify_coeffs/verify_coeff_c.f90
+!
+!! Part of the fd-schemes test suite.
+!
+! SPDX-License-Identifier: BSD-3-Clause
+
+program verify_coeff_c
+  !! Tests the computation of finite difference coefficient for non-uniform grids acting on 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 :: c
+  real, parameter :: cref = 1.0 / 9.0
+
+  real :: numerator, denominator, divisor
+  real, parameter :: numerator_f2ref = 2.0 / 3.0
+  real, parameter :: denominator_f2ref = 6.0
+  real, parameter :: divisor_f2ref = 4.0
+  
+  call initialise_suite("Verify coefficient C")
+
+  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
+
+  call coeff_c_components(points_to_deltas(stencil), numerator, denominator, divisor)
+  call test_report("Coefficient C numerator", check_scalar(numerator, numerator_f2ref * (h**3)))
+  call test_report("Coefficient C denominator", check_scalar(denominator, denominator_f2ref * (h**3)))
+  call test_report("Coefficient C divisor", check_scalar(divisor, divisor_f2ref * h))
+  c = coeff_c(stencil)
+  call test_report("Coefficient C", check_scalar(c, cref / (4.0 * h)))
+
+  call finalise_suite()
+
+end program verify_coeff_c
diff --git a/tests/fd-schemes/verify_coeffs/verify_coeff_d.f90 b/tests/fd-schemes/verify_coeffs/verify_coeff_d.f90
new file mode 100644
index 0000000..fcc22ca
--- /dev/null
+++ b/tests/fd-schemes/verify_coeffs/verify_coeff_d.f90
@@ -0,0 +1,59 @@
+! tests/fd-schemes/verify_coeffs/verify_coeff_d.f90
+!
+!! Part of the fd-schemes test suite.
+!
+! SPDX-License-Identifier: BSD-3-Clause
+
+program verify_coeff_d
+  !! Tests the computation of finite difference coefficient for non-uniform grids acting on 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 :: d
+  real, parameter :: dref = -1.0 / 9.0
+
+  real :: numerator, denominator, divisor
+  real, parameter :: numerator_f2ref =- 2.0 / 3.0
+  real, parameter :: denominator_f2ref = 6.0
+  real, parameter :: divisor_f2ref = 4.0
+  
+  call initialise_suite("Verify coefficient D")
+
+  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
+
+  call coeff_d_components(points_to_deltas(stencil), numerator, denominator, divisor)
+  call test_report("Coefficient D numerator", check_scalar(numerator, numerator_f2ref * (h**3)))
+  call test_report("Coefficient D denominator", check_scalar(denominator, denominator_f2ref * (h**3)))
+  call test_report("Coefficient D divisor", check_scalar(divisor, divisor_f2ref * h))
+  d = coeff_d(stencil)
+  call test_report("Coefficient D", check_scalar(d, dref / (4.0 * h)))
+
+  call finalise_suite()
+
+end program verify_coeff_d
diff --git a/tests/fd-schemes/verify_coeffs/verify_coeff_e.f90 b/tests/fd-schemes/verify_coeffs/verify_coeff_e.f90
new file mode 100644
index 0000000..9f5f75a
--- /dev/null
+++ b/tests/fd-schemes/verify_coeffs/verify_coeff_e.f90
@@ -0,0 +1,50 @@
+! tests/fd-schemes/verify_coeffs/verify_coeff_e.f90
+!
+!! Part of the fd-schemes test suite.
+!
+! SPDX-License-Identifier: BSD-3-Clause
+
+program verify_coeff_e
+  !! Tests the computation of finite difference coefficient for non-uniform grids acting on f_{i}.
+
+  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 :: e
+  real, parameter :: eref = 0.0
+  
+  call initialise_suite("Verify coefficient E")
+
+  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
+
+  e = coeff_e(stencil)
+  call test_report("Coefficient E", check_scalar(e, eref / (4.0 * h)))
+
+  call finalise_suite()
+
+end program verify_coeff_e
diff --git a/tests/nucfd_tests_mod.f90 b/tests/nucfd_tests_mod.f90
index fcb078d..225dfad 100644
--- a/tests/nucfd_tests_mod.f90
+++ b/tests/nucfd_tests_mod.f90
@@ -1,8 +1,12 @@
+! tests/nucfd_tests_mod.f90
+!
+!! Defines the nucfd_tests module.
+!
+! SPDX-License-Identifier: BSD-3-Clause
+
 module nucfd_tests
   !! NuCFD test module, enables defining test suites, reporting results of tests and overall status
   !! of the suite, and provides utilities for checking floating-point numbers.
-  !!
-  !! SPDX-License-Identifier: BSD-3-Clause
 
   implicit none
 
@@ -10,7 +14,8 @@ module nucfd_tests
   public :: initialise_suite, finalise_suite
   public :: test_report
   public :: check_rms
-
+  public :: check_scalar
+  
   character(len=:), allocatable :: suite_name
   logical, save :: passing
 
@@ -67,11 +72,34 @@ subroutine test_report(test_name, test_status)
 
   end subroutine test_report
 
+  logical function check_scalar(test, ref)
+    !! Compute error and report for scalar values.
+
+    real, intent(in) :: test !! The test data
+    real, intent(in) :: ref  !! The reference data
+
+    real :: err
+    logical :: test_passing
+    
+    err = abs(test - ref)
+    if (err > (2 * epsilon(ref))) then
+       test_passing = .false.
+
+       print *, "Delta = ", err, " exceeds tolerance: ", 2 * epsilon(ref)
+       print *, "Value = ", test, " expected: ", ref
+    else
+       test_passing = .true.
+    end if
+
+    check_scalar = test_passing
+
+  end function check_scalar
+  
   logical function check_rms(test, ref)
     !! Compute RMS of error and report errors.
 
-    real, dimension(:), intent(in) :: test ! The test data
-    real, dimension(:), intent(in) :: ref  ! The reference data
+    real, dimension(:), intent(in) :: test !! The test data
+    real, dimension(:), intent(in) :: ref  !! The reference data
 
     integer :: n
     integer :: i
diff --git a/tests/tridsolver/CMakeLists.txt b/tests/tridsolver/CMakeLists.txt
new file mode 100644
index 0000000..56cb4b5
--- /dev/null
+++ b/tests/tridsolver/CMakeLists.txt
@@ -0,0 +1,23 @@
+## tests/tridsolver/CMakeLists.txt
+#
+## Description
+#
+# CMakeLists file for the tri-diagonal solver test suite.
+#
+## LICENSE
+#
+# SPDX-License-Identifier: BSD-3-Clause
+#
+
+add_library(trid_test_utils 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)
+
+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)
diff --git a/tests/tridsolver/system_00_anti-symm.f90 b/tests/tridsolver/system_00_anti-symm.f90
index 82d87d8..1b552dd 100644
--- a/tests/tridsolver/system_00_anti-symm.f90
+++ b/tests/tridsolver/system_00_anti-symm.f90
@@ -1,10 +1,14 @@
+! tests/tridsolver/system_00_anti-symm.f90
+!
+!! Defines the test for solving tridiagonal systems of anti-symmetric functions on periodic domains.
+!!
+!! Part of the tridsolver test suite.
+!
+! SPDX-License-Identifier: BSD-3-Clause
+
 program test_system_00_anti_symm
   !! Tests the solution of tridiagonal systems arising from compact finite differences of
   !! anti-symmetric functions on periodic domains.
-  !!
-  !! Part of the tridsolver test suite.
-  !!
-  !! SPDX-License-Identifier: BSD-3-Clause
 
   use nucfd_coeffs
   use nucfd_trid_solver
diff --git a/tests/tridsolver/system_00_symm.f90 b/tests/tridsolver/system_00_symm.f90
index 2bb9d8a..ce536b0 100644
--- a/tests/tridsolver/system_00_symm.f90
+++ b/tests/tridsolver/system_00_symm.f90
@@ -1,10 +1,14 @@
+! tests/tridsolver/system_00_symm.f90
+!
+!! Defines the test for solving tridiagonal systems of symmetric functions on periodic domains.
+!!
+!! Part of the tridsolver test suite.
+!
+! SPDX-License-Identifier: BSD-3-Clause
+
 program test_system_00_symm
   !! Tests the solution of tridiagonal systems arising from compact finite differences of symmetric
   !! functions on periodic domains.
-  !!
-  !! Part of the tridsolver test suite.
-  !!
-  !! SPDX-License-Identifier: BSD-3-Clause
 
   use nucfd_coeffs
   use nucfd_trid_solver
diff --git a/tests/tridsolver/system_11_anti-symm.f90 b/tests/tridsolver/system_11_anti-symm.f90
index 9ca35bf..641a35c 100644
--- a/tests/tridsolver/system_11_anti-symm.f90
+++ b/tests/tridsolver/system_11_anti-symm.f90
@@ -1,10 +1,14 @@
+! tests/tridsolver/system_11_anti-symm.f90
+!
+!! Defines the test for solving tridiagonal systems of anti-symmetric functions.
+!!
+!! Part of the tridsolver test suite.
+!
+! SPDX-License-Identifier: BSD-3-Clause
+
 program test_system_11_anti_symm
   !! Tests the solution of tridiagonal systems arising from compact finite differences of
   !! anti-symmetric functions.
-  !!
-  !! Part of the tridsolver test suite.
-  !!
-  !! SPDX-License-Identifier: BSD-3-Clause
 
   use nucfd_coeffs
   use nucfd_trid_solver
diff --git a/tests/tridsolver/system_11_symm.f90 b/tests/tridsolver/system_11_symm.f90
index d2cbcd5..fabf48f 100644
--- a/tests/tridsolver/system_11_symm.f90
+++ b/tests/tridsolver/system_11_symm.f90
@@ -1,10 +1,14 @@
+! tests/tridsolver/system_11_symm.f90
+!
+!! Defines the test for solving tridiagonal systems of symmetric functions.
+!!
+!! Part of the tridsolver test suite.
+!
+! SPDX-License-Identifier: BSD-3-Clause
+
 program test_system_11_symm
   !! Tests the solution of tridiagonal systems arising from compact finite differences of symmetric
   !! functions.
-  !!
-  !! Part of the tridsolver test suite.
-  !!
-  !! SPDX-License-Identifier: BSD-3-Clause
 
   use nucfd_coeffs
   use nucfd_trid_solver
diff --git a/tests/tridsolver/tridsol_test_utils_mod.f90 b/tests/tridsolver/tridsol_test_utils_mod.f90
index 29f6fe0..2979647 100644
--- a/tests/tridsolver/tridsol_test_utils_mod.f90
+++ b/tests/tridsolver/tridsol_test_utils_mod.f90
@@ -1,9 +1,13 @@
+! tests/tridsolver/tridsol_test_utils_mod.f90
+!
+!! Defines the tridsol_test_utils module.
+!!
+!! Part of the tridsolver test suite.
+!
+! SPDX-License-Identifier: BSD-3-Clause
+
 module tridsol_test_utils
   !! Module defining utility functions for the tridiagonal solver test suite.
-  !!
-  !! Part of the tridsolver test suite.
-  !!
-  !! SPDX-License-Identifier: BSD-3-Clause
 
   implicit none