So you have a transfer function G(s) = 1/(s^3+1) and you want to find gains P,I,D for your controller K(s) = P+D*s+I/s. As commonly known, the closed loop transfer function T(s) can be written as;
T(s) = K(s)*G(s)/(1+G(s)*K(s))
If you plugin G(s) and K(s) functions, you will obtain;
T(s) = (D*s^2 + P*s + I) / (s^4 + D*s^2 + (P+1)*s + I)
From here, it should be clear that a simple controller, K(s) = P+D*s+I/s, cannot stabilize the system (hint: routh-hurwitz stability criterion). That's why your controller needs to have a term that will show up as the s^3 term.
One option is to add a C*s^2 term to your controller. This way your K(s) = P+D*s+I/s+C*s^2, and resulting closed loop transfer function is;
T(s) = (C*s^3 + D*s^2 + P*s + I)/( s^4 + C*s^3 + D*s^2 + (P+1)*s + I )
For instance, imagine you want to place your poles to (-4+0i) (continous time model), then your characteristic equation needs to be;
(s+4)^4 = s^4 + 16*s^3 + +96*s^2 + 256*s + 256 = ( s^4 + C*s^3 + D*s^2 + (P+1)*s + I )
From here, by matching the coefficients, you can find the controller
K(s) = 255+96*s+256/s+16*s^2;
This controller stabilizes G(s), and also places all 4 of the closed loop poles to (-4+0i).
Another method to tune the gains of your controller is to perform Ruth-Hurwitz stability analysis and determine the range of controller gains that achieves closed loop stability.
These calculation are done for continous time however they are valid for discrete time as well, with some modification due to conversion from s to z domain.
Hope this helps.
s = tf('s');
G = 1/(s^3+1);
K = 255+96*s+256/s+16*s^2;
T = minreal(G*K/(1+G*K));
step(T)
