Hi Robert,
Your code is just fine except for one tiny but a curicial point that is internally built-in solver step size related issue of the control toolbox. In fact, there is only one time delay. Here is the correct code that addresses solver step size in computing the systems' step responses. Now, your simulation results of the both systems converge perfectly well.
clc; clear variables
%
C = tf([1.5, 47.12], [1.0, 0.0]);
G_delay = tf(1, 1, 'OutputDelay', 1.0 / 40.0e3);
P = tf([6773, 5.131e7], [1.0, 1.123e4, 5.193e7]);
H = tf(5.988e19, [1.0, 348900, 4.591e10, 2.7e15, 5.988e19]);
% build closed-loop using feedback()
T_CL = minreal(feedback(P * C * G_delay, H));
% build closed-loop manually
G_y_r = minreal((P * C * G_delay) / (1.0 + P * C * G_delay * H));
%
disp('System T_CL: '); T_CL %#ok
disp('System G_y_r: '); G_y_r %#ok
%
figure;
hold on;
dt = 1/1000; % Time step
t=0:dt:.5; % Simulation time space
step(T_CL, t, '-b');
step(G_y_r, t, '*r');
grid on;
legend('T_{CL}', 'G_{y, r}', 'Location', 'best');
Good luck.
