diff --git a/pymatgen/phonon/dos.py b/pymatgen/phonon/dos.py
index dbe5d579eea..0c06a35393e 100644
--- a/pymatgen/phonon/dos.py
+++ b/pymatgen/phonon/dos.py
@@ -354,6 +354,47 @@ def mae(self, other: PhononDos, two_sided: bool = True) -> float:
 
         return self_mae
 
+    def get_last_peak(self, threshold: float = 0.1) -> float:
+        """Find the last peak in the phonon DOS defined as the highest frequency with a DOS
+        value at least threshold * height of the overall highest DOS peak.
+        A peak is any local maximum of the DOS as a function of frequency.
+        Use dos.get_interpolated_value(peak_freq) to get density at peak_freq.
+
+        TODO method added by @janosh on 2023-12-18. seems to work well in most cases but
+        was not extensively tested. PRs with improvements welcome!
+
+        Args:
+            threshold (float, optional): Minimum ratio of the height of the last peak
+                to the height of the highest peak. Defaults to 0.1 = 10%. In case no peaks
+                are high enough to match, the threshold is reset to half the height of the
+                second-highest peak.
+
+        Returns:
+            float: last DOS peak frequency (in THz)
+        """
+        first_deriv = np.gradient(self.densities, self.frequencies)
+        second_deriv = np.gradient(first_deriv, self.frequencies)
+
+        maxima = (  # maxima indices of the first DOS derivative w.r.t. frequency
+            (first_deriv[:-1] > 0) & (first_deriv[1:] < 0) & (second_deriv[:-1] < 0)
+        )
+        # get mean of the two nearest frequencies around the maximum as better estimate
+        maxima_freqs = (self.frequencies[:-1][maxima] + self.frequencies[1:][maxima]) / 2
+
+        # filter maxima based on the threshold
+        max_dos = max(self.densities)
+        threshold = threshold * max_dos
+        filtered_maxima_freqs = maxima_freqs[self.densities[:-1][maxima] >= threshold]
+
+        if len(filtered_maxima_freqs) == 0:
+            # if no maxima reach the threshold (i.e. 1 super high peak and all other peaks
+            # tiny), use half the height of second highest peak as threshold
+            second_highest_peak = sorted(self.densities)[-2]
+            threshold = second_highest_peak / 2
+            filtered_maxima_freqs = maxima_freqs[self.densities[:-1][maxima] >= threshold]
+
+        return max(filtered_maxima_freqs)
+
 
 class CompletePhononDos(PhononDos):
     """This wrapper class defines a total dos, and also provides a list of PDos.