A quick look at this shows a common factor of 10 between the first low-frequency pole (1.3 rad/s) and the first zero at 13 rads/s (slope is now 0) then the second zero (slope is +1) with a final second pole to set the slope back to 0. The below curves seem ok then:
You can see the first pole when the phase crosses -45°, then the zero is not far otherwise the phase would further drop.
The question is how did I get the first pole and the zero? I read the curve starting from the flat portion (0 dB) and then tried to identify the -3-dB point on the magnitude graph or the -45° on the phase plot. It is an approximation here and implies that the pole and the zero are well spread. Then, once the pole is obtained, I determined the magnitude of a pole-zero pair. I then solved for the zero position which brings a valley in the magnitude plot at -19.7 dB roughly read from the magnitude curve. Voilà !
You wrote the transfer function the right way, by keeping the poles and zeros well factored. This is what is called a low-entropy format in which a leading term would indicate what the gain is for \$s=0\$ in your case (1 or 0 dB). Writing these transfer functions the right way is part of the fast analytical circuits techniques or FACTs introduced many years ago. They describe how to determine transfer functions in a swift and efficient manner and how to format a transfer function so that poles, zeros and gain immediately show up when looking at the equation.