Bode plot from from step response data
First of all - I'm not familiar with Python and didn't check your code.
1a) For FFTs it is the usual that the frequency resolution equals the inverse length of your data record.
1b) As far as I understand your code (I may be wrong here), you take the FFT of the step response. The Bode diagram usually shows the FFT of the impulse response. Take care - the result may not be what you expect.
1c) In practice (due to transients, noise,...) the usual approach to measure the FRF is

excite the system with a broadband signal, estimate the spectral density of the input and the output signal, divide the output spectral density by the input spectral density. This approach is well covered in the literature (and called &quot;spectral method&quot; or something similar

exite the system with a sinuosoid, wait until the transient has died out, and read the magnitude and phase from your data. Repeat for every frequency where you want to know the FRF.

There are other approaches to directly measure the FRF, like the Local Polynomial Method (or Local Rational Method), TRIMM, and more - they are covered in (scientific) publications, together with a lot of beautiful maths. The approaches I mentioned are called non-parametric methods (there are also parametric methods, but that may lead to far away for your use case)



The phase jump at Nyquist frequency is to be expected. The FFT spectrum of a rational signal will be symmetric to the DC-frequency - the magnite is axisymmetric, the phase point symmetric. For some FFTs the &quot;negative&quot; and &quot;positive&quot; frequencies are wrapped (this may depend on your FFT implementation). The result is a &quot;jump&quot; the phase.

The input signal should excite the entire frequency range of interest - check the spectral density function of the input signal (i.e. the step).