it sounds like there are two issues here: 1. estimating a frequency response using test inputs, 2. what does lowpass do. Let's tackle them in sequence.
Estimating the frequency response: You were on the right track with injecting sine waves of different frequencies and comparing the input to the output. A couple of things to keep in mind. You need to make sure to either delete the initial transient from the output (and the corresponding input time window in the input) or ensure that the initial transient is a small proportion of the overall response. Keep in mind that the response you're trying to estimate is for a particular frequency. Taking the ratio of the RMS's to estimate the magnitude response might be o.k. for linear systems, but in general you have to be careful that you're not RMS'ing components of the output at frequencies other than the input test frequency.
To get the phase response in the time domain you need to estimate the delay between the input and the output at the test frequency and then convert to phase. For sure, taking the angle of the RMS ratio will not yield that delay. If you really want to estimate the delay in the time domain, there is a function called finddelay that may be of use.
However, the alternative approach is to work in the frequncy domain, i.e., use the ratio of the frequency response of the output to the frequency resposne of the input to directly stimate the magnitude and phase response of the system. The code that follows shows how to do that for a simple low pass filter. Keep in mind, that any approach you use may begin to suffere as you test frequency gets close to the Nyquist frequency. You can experiment with this code to see how close you can get before this simple approach begins to break down.
fs = 250;
Ts = 1/fs;
[b,a] = butter(3,0.25);
fplot = logspace(-2,pi,1000)/pi*fs/2;
hplot = freqz(b,a,fplot,fs);
f = 30;
t = 0:Ts:(100/f); t = t(:);
u = sin(2*pi*f*t);
y = filter(b,a,u);
f_fft = (0:length(t)-1)/(length(t)-1)*fs;
Y_fft = interp1(f_fft,fft(y),f);
U_fft = interp1(f_fft,fft(u),f);
h = freqz(b,a,[0 f],fs); h = h(2);
m_comp = 20*log10([abs(h) abs(Y_fft/U_fft) (rms(y)/rms(u))])
p_comp = 180/pi*angle([h Y_fft/U_fft])
As you discovered, the low pass filter returned from the lowpass function you called is indeed a linear phase filter. However, the lowpass function applies that filter in such a way to yield something very close to zero-phase filtering, which is what I think the doc means by "compensates for the delay introduced by the filter," which is not a very clear statement and is not explained in the doc. The filtfilt function does something similar, and is explained in doc filtfilt. Let's compare the output of lowpass for a particular input to the output of filtfilt using the filter returned by lowpass. Note that the effective gain of filtfilt is the magnitude squared of the filter, so the output of filtfilt has to be divided by the filter gain to match the result from lowpass. If you zoom in on the plot, you'll see that lowpass and filtfilt must use different approaches near the intial and final times of the response.
t = 0:1/fs:10-1/fs;
x = sin(2*pi*t'*f);
y = filtfilt(d.Coefficients,1,x(:,end));
h = freqz(d,[0 f],fs); h = h(2);