How to do the bode plot for this function without Matlab?
When you have such transfer function (TF), it is good, in my opinion, to rewrite it in a low-entropy format in which poles and zeroes explicitly appear. For instance, you can see that this TF is made of a second-order TF cascaded with that of a zero and a pole located at the same position. The pole and the zero cancel each others in magnitude. Therefore, the magnitude response is that of a second-order TF: flat then going down with a -2-slope. We will see the phase later on.
Let's reformat this equation:

From here, you see a dc gain (actually an attenuation), then a normalized second-order polynomial and, finally, the pole and the zero. The negative sign in the zero indicates that it is in the right half-plane: this is a RHP zero whose magnitude is that of a classical zero but it lags the phase (like a pole would) rather than leading it. So the goal is to factor the whole thing as follows:

If you do the maths ok, you should find:

So you have a straight line starting at -23.5 dB in the magnitude plot, flat until 3.9 rd/s, then a slight peak occurs (\$Q\$ is almost 1) followed by a -2-slope (-40 dB per decade). As I said, the pole and the RHP zero are both coincident and neutralize each other in magnitude. The phase will start from 180° (you have an inversion) and will asymptotically hit 0°. The pace at which the phase drops depends on \$Q\$ and it should be rather smooth with 1. Then, the pole and zero phase will now cumulate (remember, this is a RHP zero) and prolong the phase drop from 0° to -180°. The plot could look like this:

Without the pole-zero pair, the phase of the inverting second-order network would cross 90° at the resonance and further go down to 0°. Adding the pole-zero pair positioned before the resonance will distort this value and I have indicated an approximate point (it's 3.08 rd/s) at which crossing occurs but for a hand-drawn response it should be good enough. Finally, using Mathcad, you obtain the exact shape: