Performing sine sweep for system identification

I am working on a project in which I am provided with a plant that is given an input over a period of time and creates a vector of outputs for each instance in time that is fed into the plant. My goal is to perform sine sweep with the input and outputs of the plant to fit a transfer function to the plant. I cannot go into the plant in Simulink as I am provided with only executable P-code. I am having trouble running my code for a number of frequencies as it is taking a very long time, and creating a bode plot that does not look right because I am checking it using the iddata command with the input and output vectors from the plant with the respective sampling period. I have attached the files necessary to run the plant in a zipped folder as well as included my code below. If anyone could offer suggestions on improving my code or approaching the sine sweep method differently it would be greatly appreciated. Thank you.

clc;
clear;
N=25;
w=logspace(-1,3,N);
g=zeros(length(w),1);
gdb=zeros(length(w),1);
phi=zeros(length(w),1);
T=zeros(length(w),1);
A=1;

for ii=1:length(w)
    
    T(ii)=(2*pi)/(w(ii));
%     T(ii)=P(ii);
    modelName = 'openLoopTestBed_R2019b';
    Ts=0.01; %sampling period, in seconds (rather, time units)
    tfinal = 100;
    timeInput = [0:Ts:tfinal]';
    uValues = A*sin(w(ii)*timeInput);
    u = [];
    u.time = timeInput;
    u.signals.values = uValues;
    u.signals.dimensions = 1;
    Out = sim(modelName,'StopTime',num2str(tfinal));
    timeOut = Out.y.time;
    yValues = Out.y.signals.values;
        
    Y_inphase=yValues.*sin(w(ii)*timeInput);
    Y_quadrature=yValues.*cos(w(ii)*timeInput);
    Ip=(1/T(ii))*Ts*trapz(Y_inphase);
    Iq=(1/T(ii))*Ts*trapz(Y_quadrature);
    g(ii)=(2/A)*sqrt(Ip^2+Iq^2);%gain
    phi(ii)=atan2(Iq,Ip);%phase
    
end
gdb=20*log10(g);
phi=unwrap(phi)*(180/pi);
subplot(2,1,1)
semilogx(w,gdb)
subplot(2,1,2)
semilogx(w,phi)

idDataForThePlant = iddata(yValues,uValues,Ts);
Phat_DT = etfe(idDataForThePlant);
figure
bode(Phat_DT)

I cannot look at the information for the plant

[Pat Gipper]

The plant response has sizable low-frequency noise added to it that is going to wreck havoc on a linear analyzer unless you run a lot of cycles for each frequency input. Just running a series of step responses for about 500 seconds allowed me to "average" out the noise and I see a linear gain of 3 and approximate time constant of 20 seconds. Also it seems to saturate with inputs larger than +/-10.