Simulated plots of the reflection transfer function (magnitude and phase) of a two-mirror Fabry-Perot cavity, as well as a simulated PDH signal.
For the cavity mirrors I chose amplitude reflectivities of r1 = 0.99 and r2 = 0.98, and a length L = 1 m. For the modulation frequency I chose fm = 23 MHz (not 25 MHz as written in upload log).
Matplotlib code is as follows:

from matplotlib import plt
import numpy as np
c = 3e8
L = 1.
fsr = c/(2*L)
f = np.arange(-.4*fsr,.4*fsr,fsr/1e4)
r1 = .99
r2 = .98
t1 = (1-r1**2)**.5
R = (r1-(r1**2+t1**2)*r2*np.e**(2j*2*np.pi*f*L/c))/(1-r1*r2*np.e**(2j*2*np.pi*f*L/c))
fm = 23e6
Rffm = (r1-(r1**2+t1**2)*r2*np.e**(2j*2*np.pi*(f+fm)*L/c))/(1-r1*r2*np.e**(2j*2*np.pi*(f+fm)*L/c))
Rfnfm = (r1-(r1**2+t1**2)*r2*np.e**(2j*2*np.pi*(f-fm)*L/c))/(1-r1*r2*np.e**(2j*2*np.pi*(f-fm)*L/c))
pdh = R*np.conjugate(Rffm)-np.conjugate(R)*Rfnfm
plt.subplot(311)
plt.plot(f/1e6,100*np.abs(R)**2,'#880088')
plt.axis([-.4*fsr/1e6,.4*fsr/1e6,0,102])
plt.ylabel('Reflected power (\%)')
plt.subplot(312)
plt.plot(f/1e6,180*np.angle(R)/np.pi,'#880088')
plt.axis([-.4*fsr/1e6,.4*fsr/1e6,-35,35])
plt.ylabel('Reflected phase (deg.)')
plt.subplot(313)
plt.plot(f/1e6,np.imag(pdh),'#880088')
plt.ylabel('PDH readout (arb.)')
plt.xlabel('$f-f_\mathrm{res}$ (MHz)')
This  plot was created with Matplotlib.
Zueignung