syms z;
H = (1-z^-1) / (1-3*z^-1+2*z^-2);
[N, D] = numden(H);
Nc = eval(coeffs(N)); %Get coeffs and evaluatle symbolic variable, i.e. make real matrix
Dc = eval(coeffs(D));
Nc = Nc./(Dc(1)); %Turn into proper polynomial, first coeffs of a is 1
Dc = Dc./(Dc(1));
M = idpoly(Dc, Nc, 'NoiseVariance',0)