diff --git a/CHANGELOG.md b/CHANGELOG.md
index 53d7baa9..1628adee 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -50,6 +50,8 @@ This project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.htm
 
 - Running the configure command now requires the correct combination of library locations to be specified; \[Brandon Lazard; 6-21-2024] [#535](https://github.com/geodynamics/Rayleigh/pull/535) 
 
+- Fixed a bug that was causing Rayleigh to crash during output if finite-difference mode was active and second-order derivatives in radius were required for output. \[Nick Featherstone; 6-26-2024; [#566](https://github.com/geodynamics/Rayleigh/pull/566) \]
+
 
 ## [1.2.0] - 5-29-2024
 
diff --git a/src/Diagnostics/Diagnostics_Second_Derivatives.F90 b/src/Diagnostics/Diagnostics_Second_Derivatives.F90
index 62fa70b9..72f9519f 100644
--- a/src/Diagnostics/Diagnostics_Second_Derivatives.F90
+++ b/src/Diagnostics/Diagnostics_Second_Derivatives.F90
@@ -27,6 +27,7 @@ Module Diagnostics_Second_Derivatives
     Use Diagnostics_Base
     Use Structures
     Use Spectral_Derivatives
+    Use Finite_Difference, Only : d_by_dx3d3
     Implicit None
 
 
@@ -378,8 +379,12 @@ Subroutine Compute_Second_Derivatives(inbuffer)
         !Move to p1b configuration
         Call d2buffer%reform()
         Call d2buffer%construct('p1a')
-        Call gridcp%To_Spectral(d2buffer%p1b,d2buffer%p1a)
-        Call gridcp%dealias_buffer(d2buffer%p1a)
+        If (chebyshev) Then
+            Call gridcp%To_Spectral(d2buffer%p1b,d2buffer%p1a)
+            Call gridcp%dealias_buffer(d2buffer%p1a)
+        Else
+            d2buffer%p1a = d2buffer%p1b
+        Endif
         d2buffer%p1b = 0.0
         d2buffer%config='p1a'
 
@@ -387,11 +392,18 @@ Subroutine Compute_Second_Derivatives(inbuffer)
 
         !////////////////////////////////////////////////////////////
         ! Steps 5-6:  Load d2_by_dr2 and d2_by_drdt into the buffer
-        Do i = 1, nddfields*2
-            j = i+nddfields*2
-            Call gridcp%d_by_dr_cp(i,j,d2buffer%p1a,1)
-        Enddo
-        Call gridcp%From_Spectral(d2buffer%p1a,d2buffer%p1b)
+        If (chebyshev) Then
+            Do i = 1, nddfields*2
+                j = i+nddfields*2
+                Call gridcp%d_by_dr_cp(i,j,d2buffer%p1a,1)
+            Enddo
+            Call gridcp%From_Spectral(d2buffer%p1a,d2buffer%p1b)
+        Else
+            Do i = 1, nddfields*2
+                j = i+nddfields*2
+                Call d_by_dx3d3(i,j,d2buffer%p1a,1)
+            Enddo
+        Endif
         d2buffer%p1a=d2buffer%p1b
         Call d2buffer%deconstruct('p1b')
 
diff --git a/src/Makefile.fdeps b/src/Makefile.fdeps
index 0b7dbf82..6530efb2 100644
--- a/src/Makefile.fdeps
+++ b/src/Makefile.fdeps
@@ -25,7 +25,7 @@ Diagnostics_Mean_Correction.o : Diagnostics_Mean_Correction.F90 indices.F Diagno
 Diagnostics_Miscellaneous.o : Diagnostics_Miscellaneous.F90 indices.F Diagnostics_Base.o 
 Diagnostics_Poynting_Flux.o : Diagnostics_Poynting_Flux.F90 indices.F Diagnostics_Base.o 
 Diagnostics_Scalars.o : Diagnostics_Scalars.F90 indices.F Controls.o Diagnostics_Base.o 
-Diagnostics_Second_Derivatives.o : Diagnostics_Second_Derivatives.F90 indices.F Spectral_Derivatives.o Structures.o Diagnostics_Base.o 
+Diagnostics_Second_Derivatives.o : Diagnostics_Second_Derivatives.F90 indices.F Finite_Difference.o Spectral_Derivatives.o Structures.o Diagnostics_Base.o 
 Diagnostics_Thermal_Energies.o : Diagnostics_Thermal_Energies.F90 indices.F Diagnostics_ADotGradB.o Diagnostics_Base.o 
 Diagnostics_Thermal_Equation.o : Diagnostics_Thermal_Equation.F90 indices.F Diagnostics_ADotGradB.o Diagnostics_Base.o 
 Diagnostics_Thermodynamic_Gradients.o : Diagnostics_Thermodynamic_Gradients.F90 indices.F Diagnostics_Base.o