How to calculate poles for a high order system?
A second order system is relatively simple. It is straight forward to determine the overshoot and the settling time. This is not the case for a higher order system.

I will limit the discussion to linear systems having linear controllers.

For a higher order system you generally construct a cost function. There are many ways to do so. A good place to start is with the linear quadratic regulator (lqr in MATLAB). For your SISO system it will have the form
$$\int_0^\infty (x^\top Q x + R u^2 + 2 x^\top N u) \ \mathrm{d} t.$$
A good place to start is to set
$$\int_0^\infty (x^\top x + R u^2) \ \mathrm{d} t.$$
Then you can vary R until you get a satisfactory response. The lqr function in MATLAB will give you the feedback matrix.

Because you are after low overshoot and fast settling times the LQR is not really the best tool you can use. Instead the ITAE (integrated time absolute error) minimizes
$$\int_0^\infty |e| t \ \mathrm{d}t.$$
This way you penalize errors more the further away they occur.

For a fourth order system your target transfer function is
$$\frac{\omega_0^4}{s^4 + 2.1 \omega_0 s^3 + 3.4 \omega_0^2 s^2 + 2.7 \omega _0^3 s + \omega_0^4}.$$
You can find more information on page 21 here.

Setting omega to 10 yields


For your system you do not need to control the zeros because the s term is small. If it were larger, you would need to remove the zeros.