Nice function especially useful for reducing systems coming from matlab linearizations.
There is a little mistake in the code easily catch by M-lint
wrong
if ~isreal(r) || (fix(r)~=r) || (r<1) || (r>n) %%% here n is undefined
error('Invalid value of reduced model order')
end
[num,den]=tfdata(G,'v');
D_fact=num(1)/den(1);
num=num-D_fact*den;
num1=num(end:-1:1)/den(1);
den1=den(end:-1:1)/den(1);
n=length(den1)-1;