PI Controller Design (Ziegler-Nichols) and Unit Step Response Graphing for 4th order system
From what I can see the recipe table is not telling the whole truth about Kp. It is missing gain of the step response. This expression is valid as long as the step response is tending to the same value as the input step therefore their gain ratio is 1. But system G is responding to a unit step by rising to 2. So step/response ratio of the given system is 1/2.

Kp for PI controllers with Ziegler-Nichols First Method in fact is
\begin{equation}
K_p=0.9⋅\frac{K_{step}}{K_0}⋅\frac{T}{L}
\end{equation}
Because gain of a unit step is 1 this is
\begin{equation}
K_p=0.9⋅\frac{1}{K_0}⋅\frac{T}{L}
\end{equation}
Only with the step response tending to 1 so its overall gain K0 is equal to 1 as well this eventually becomes
\begin{equation}
K_p=0.9⋅\frac{1}{1}⋅\frac{T}{L}=0.9\frac{T}{L}
\end{equation}
as seen in the recipe table. But that is not the case with given sytem G. So gain ratio has to be considered as well as time ratio when obtaining Kp of G(t):
\begin{equation}
K_p=0.9⋅\frac{1}{K_0}⋅\frac{T}{L}=0.9⋅\frac{1}{2}⋅\frac{2.5188}{0.441}\approx2.57
\end{equation}
(Source: https://onlinelibrary.wiley.com/doi/pdf/10.1002/9780470029558.oth5 , page 292)
Also the penultimate line of the Matlab code
Gs=feedback(Gr*Css,1);

is looping 'Gr' which is the given system G(s) and 'Css' which is the FOPDT approximant for this very system G(s). Looping G(s) and the PI-controller 'Cpi' instead should to the trick.
Gs=feedback(Gr*Cpi,1);

Still that closed loop is unstable as it has two positive conjugated complex poles:
pole(Gs)
ans =
  -4.0152 + 2.0976i
  -4.0152 - 2.0976i
   0.3426 + 2.0134i   % uh
   0.3426 - 2.0134i   % oh
  -0.6548 + 0.0000i

Using Photoshop for the graphical approximation (pixel perfect mesurement because why not) parameters of the FOPDT approximant are more like:
$$K_0=2$$
$$L=0.7$$
$$t=2.26$$
With these the step response is a stabilizing underdamped oscillation.

New Matlab code:
s=tf('s');%Time Transfer function

Gr=32/(s+2)^4;%system G(s)

% Parameters of FOPDT approximant
K0=2; %step response gain
L=.7; %step response delay time
t=2.26; %step response time constant estimate

% Ziegler-Nichols First Method
% PI controller
Ti=L/.3; %integral time
Kp=.9/K0*t/L; %proportional gain


Cpi=Kp*(1+1/(Ti*s));%transfer function of PI


%closed loop of G(s) and PI
Gs=feedback(Gr*Cpi,1);

step(Gs)