Fp = [4e3 8e3];
Wp = 2*pi*Fp;
Fs = [1e3 16e3];
Ws = 2*pi*Fs;
Rp = 1;
Rs = 20;

[N,Wn]=buttord(Wp,Ws,Rp,Rs,'s');
[b,a]=butter(N,Wn,'s');
subplot(221)
plot(roots(a),'x')
hold on
plot(roots(b),'o')
grid
axis([-1 1 -1 1]*1e5)

subplot(222)
F = 0:10:40e3;
H = freqs(b,a,2*pi*F);
plot(F/1000,20*log10(abs(H)));
grid
axis([0 20 -60 0])

subplot(223)
plot(F/1000,20*log10(abs(H)));
axis([0 10 -2 0])
grid

subplot(224)
plot(F/1000,20*log10(abs(H)));
axis([0 20 -30 -10])
grid