Replicating the PSD values from plot_psd

I am trying to compare the PSDs between multiple EDF files on a series of channels. First I want to calculate the values that are plotted when calling plot_psd. Then I want to find the area under the curve for certain ranges (alpha, beta, theta, gamma, etc) for comparison (should be trivial once I get the first step down). Below is the code I used for my calculations. You can see the two outputted plots look pretty different. How can I ensure I get the same values and achieve my desired values? Thanks for any help.

raw ='file.edf', exclude=ignored_channels)
spectra, freqs = psd_multitaper(raw, fmin=0.5, fmax=40, tmin=0, tmax=250)

plt.plot(freqs.tolist(), spectra[0], linewidth=0.25)

raw.plot_psd(fmin=0.5, fmax=40, tmax=250, tmin=0)

My Calculations

MNE-Python version: 0.23.0
Operating System: macOS 11.4

raw.plot_psd() uses mne.time_frequency.psd_welch(). Furthermore, you are plotting the PSD with dB=True, which you are not doing in your manual calculation.