Need for phase unwrapping.

[faze-geek]

There is a discrepancy between SymPy and other control toolkits concerning the phase data generated by the bode plot for frequency response.
By definition , phase = arg(H(jw)) which is correctly implemented. The difference in value arrives due to the the principal argument of arg being (180,180].
Due to this there is a sudden change in behavior of the phase plot as values tend to -180 (sudden jump by ~360 which is not the true nature of the system). To make the phase plot continuous, we would require phase unwrapping.

Here is an example -
Sympy

>>> tf = TransferFunction(1,s**3 + 2*s**2 + s,s)
>>> bode_phase_plot(tf) # tends to -180 degree at 1

[faze-geek]

See how other Control Toolkits offset that jump.

MATLAB

>> H = tf([1],[1 2 1 0])
>> bode(H)

python-control

G = (1)/(s**3 + 2*s**2 + s)
bode_plot(G)

[sylee957]

For nontriviality, there are fully algebraic method to compute the phase unwrapping
You may have to consult some research like in https://doi.org/10.1109/78.678480