You will need to use the System Identification Toolbox for this. If you have both the input and output time-domain signals, start wiith the iddata function and then (preferably) use the ssest function to estimate the system. Use the compare function to see how well the estimiated system matches the data.
You can get a quick estimate of the system dynamics by plotting the imaginary part of the Fourier transform of the transfer function (the Fourier transform of the output signal divided element-wise by the Fourier transform of the input signal). The poles are the frequenciies at which it approaches infinity, and the zeros are where the curve crosses zero. That should give you an estimate of the system order that you can use wiith ssest.
.
