Determine filter function from input and output signals

I have an known input signal that I am putting through a "black box" resulting in an experimentally measured output signal. I want to model the "black box" by deducing its filter function from the input and output signal. How can I do that?

The code below is a test code to develop the filter deduction capabilities. First I generate a square wave. Then I create the output signal using a first order butterworth filter. This is to simulate the "black box". Then I use invfreq to try to deduce the filter form. However, it does not recreate the orignal butterworth filter.

clear
close all

%% Generate Signal

%sampling frequency
f_sample = 200e6; %samples/s
%sampling period
tau_sample = 1/f_sample; %s

%time vector
t = 0:tau_sample:500e-6; %s

%vector length
L = length(t);
%frequencies
f = (f_sample*(0:(L/2))/L)';

%wave period
tau_wave = 100e-6; %s

%input square wave
vin = square(t*2*pi/(tau_wave));


%% Filter Input signal through "Black Box" to create output signal

%nyquist sampling frequency
Ny = f_sample/2;
%order
n=1;
%cutoff frequency
LP = 14e3;
Wn = LP/Ny;
[b1,a1] = butter(n, Wn,'low');
H1=dfilt.df2t(b1,a1);
HGfilter=H1;

%frequency response of filter
[h,w] = freqz(b1,a1,length(f),f_sample);

%output signal after filter
vout = filter(HGfilter,vin);


%% Fourier Transform

% raw FFT
Vin = fft(vin);

%normalize FFT coefficients
P2 = abs(Vin/L);
Pin = P2(1:L/2+1);
Pin(2:end-1) = 2*Pin(2:end-1);

% raw FFT
Vout = fft(vout);

%normalize FFT coefficients
P2 = abs(Vout/L);
Pout = P2(1:L/2+1);
Pout(2:end-1) = 2*Pout(2:end-1);

%% Deduced filter

%filter transmission
T = abs(Pout./Pin);
%calc transfer function coefficients
[b,a] = invfreqz(T,f,2,3);
%transfer funciton
hz = tf(b,a);
%frequency response
[h_calc, w_calc] = freqz(b,a,length(f),f_sample);
%model filtered signal
[v_calc,t_calc] = lsim(hz,vin,t);

% raw FFT
Vcalc = fft(v_calc);

%normalize FFT coefficients
P2 = abs(Vcalc/L);
Pcalc = P2(1:L/2+1);
Pcalc(2:end-1) = 2*Pcalc(2:end-1);

%%  Plot

figure(101)
set(gca, 'FontSize', 18);
plot(t, vin, 'LineWidth', 2)
hold on
plot(t, vout, 'LineWidth', 2);
plot(t_calc,v_calc, 'LineWidth', 2);
hold off
%axis([2.42 2.52 -0.7 10]);
xlabel('Time (\mus)','FontSize', 20, 'FontWeight', 'bold');
ylabel('Signal (a.u.)', 'FontSize', 20, 'FontWeight', 'bold');
legend('Input', 'Measured', 'Filtered');


figure(102)
set(gca, 'FontSize', 18);
plot(f, Pin, 'LineWidth', 2)
hold on
plot(f, Pout, 'LineWidth', 2);
plot(f, Pcalc, 'LineWidth', 2);
hold off
%axis([2.42 2.52 -0.7 10]);
xlabel('Frequency (Hz)','FontSize', 20, 'FontWeight', 'bold');
ylabel('FT (a.u.)', 'FontSize', 20, 'FontWeight', 'bold');
legend('Input', 'Measured','Filtered');
set(gca,'XScale','log')
set(gca,'YScale','log')

%plot filter response
figure(103)
plot(w,h)
hold on
plot(w_calc,h_calc)
hold off
xlabel('Frequency (Hz)','FontSize', 20, 'FontWeight', 'bold');
ylabel('Transfer Function (a.u.)', 'FontSize', 20, 'FontWeight', 'bold');
legend('True Filter','Deduced Filter')
set(gca,'XScale','log')
set(gca,'YScale','log')