function wilkinsonpol

p=poly(1:20);
q=p;
r=p;
s=p;
q(21)=p(21)*(1+10^-5);
r(11)=p(11)*(1+10^-5);
s(6)=p(6)*(1+10^-5);

rp=roots(p);
rq=roots(q);
rr=roots(r);
rs=roots(s);

hold on
h1=plot(real(rp),imag(rp),'b.','Markersize',15);
h2=plot(real(rq),imag(rq),'r.','Markersize',15);
h3=plot(real(rr),imag(rr),'g.','Markersize',15);
h4=plot(real(rs),imag(rs),'k.','Markersize',15);
legend('(a) no perturbation','(b) a_1 perturbed','(c) a_{10} perturbed','(d) a_{15} perturbed');
hold off