frequency response of a time domain subroutine in matlab
The easiest would be to do an FFT of the impulse response. Something like this
%% calculating the transfer function of a time domain LTI model

% create an FFT grid
fftLength = 1024;
% make a unit impulse
diracDelta = zeros(fftLength,1); diracDelta (1) = 1;
% run your &quot;black box&quot; function, to create the impulse response
h = myFunction (diracDelta); % plus whatever arguments you need
% FFT
transferFunction = fft(h);

You can certainly encapsulate this in a function but you need to pass in your black box function as a function handle and need to figure out how to manage additional arguments and state keeping (for filters with state). Given that the code above is so simple, it doesn't seem worth the bother.
Keep in mind that this ONLY works for LTI system. For anything else, the transfer function is undefined.
The choice of the FFT length is important: to avoid time domain aliasing it needs to be long enough to contain the entire impulse response or at least where it's decayed to &quot;low enough&quot;.