Hi
Your basic procedure is correct, but you mixed up the parameters a bit. I admit that this is a bit confusing, but the structure that the Ziegler-Nichols table is written for is
C(s) = Kp*(1+1/(s*Ti)+s*Td)
I.e. you get Ti and Td, and not Ki and Kd.
Your measured values seem a bit off btw, did you read them from the plot?
You can check the exact values with e.g. the margin command, which gives you the gain margin and the corresponding frequency (i.e. your Ku and Tu). Note that the frequency will be given in rad/s and not Hz, so you need to convert that.
You should also be aware that ZN is a heuristic method, and will not work or return meaningful results for arbitrary plants.
Here a small script to compute the Ziegler Nichols Parameters for a system (programmatically, not with simulink):
clear all
close all
clc
M = 1;
D = 8;
K = 100;
s = tf('s');
P = 1/(s+1)^5;
[ Ku, ~, WKu, ~ ] = margin(P);
fKu = WKu/(2*pi);
Tu = 1/fKu;
fprintf('Limit Kp: %g, Period: %g [s]\n',Ku,Tu)
Kn = 0.01;
Kp0 = 0.4;
Ki0 = 0.1;
Kd0 = 0.1;
C_PID_0 = pid(Kp0,Ki0,Kd0,Kn);
Kplim = 0.9*Ku;
Kilim = 0;
Kdlim = 0;
C_PID_lim = pid(Kplim,Kilim,Kdlim,Kn);
KpZN = 0.6*Ku;
TiZN = 0.5*Tu;
KiZN = KpZN/TiZN;
TdZN = 0.125*Tu;
KdZN = KpZN*TdZN;
C_PID_ZN = pid(KpZN,KiZN,KdZN,Kn);
fprintf('Ziegler Nichols Parameters:\n')
fprintf('Kp: %g, Ki: %g [s], Kd: %g [s]\n',KpZN,KiZN,KdZN);
G0 = feedback(P*C_PID_0,1);
G1 = feedback(P*Kplim,1);
G2 = feedback(P*C_PID_ZN,1);
figure(1)
nyquist(P*C_PID_0,P*C_PID_lim,P*C_PID_ZN)
grid on
axis([ -2.8419 3.6822 -4.7619 4.3024 ])
legend('low gains','limit','Ziegler Nichols')
figure(2)
step(G0,G2)
ylim([- 2 2 ])
legend('Low Gains','Ziegler Nichols','Location','SouthEast')
Cheers
EDIT: typo in ZN gain calculation, fixed