  
close all
format compact
n = 20
syms x
WP = prod(x-(1:n)); %Konstruerea Wilkinsons polynom
% WP = 1*x^20 - 210*x^19 + 20615*x^18 - 1256850*x^17 + 53327946*x^16 - 1672280820*x^15 + 40171771630*x^14 - 756111184500*x^13 + 11310276995381*x^12 -...
%      135585182899530*x^11 + 1307535010540395*x^10 - 10142299865511450*x^9 + 63030812099294896*x^8 - 311333643161390656*x^7 + 1206647803780373248*x^6 - ...
%      3599979517947607040*x^5 + 8037811822645051392*x^4 - 12870931245150988288*x^3 + 13803759753640704000*x^2 - 8752948036761600000*x + 2432902008176640000;
Z = sort(solve(WP)); %Hitta rötterna till Wilkinsons polynom exakt
plot(real(Z),imag(Z),'o','linewidth',2) %Rita upp rötterna, skall vara 1 tom 20
xlabel('RE')
ylabel('IM')
hold on

for i = 1:20
WP_perturbed = WP + randn*1e-10*x^20 + ... %Stör koefficienterna till Wilkinsons polynom med ett slumpmässigt tal av storleksordning 10^(-10)
                    randn*1e-10*x^19 + ...
                    randn*1e-10*x^18 + ...
                    randn*1e-10*x^17 + ...
                    randn*1e-10*x^16 + ...
                    randn*1e-10*x^15 + ...
                    randn*1e-10*x^14 + ...
                    randn*1e-10*x^13 + ...
                    randn*1e-10*x^12 + ...
                    randn*1e-10*x^11 + ...
                    randn*1e-10*x^10 + ...
                    randn*1e-10*x^9 + ...
                    randn*1e-10*x^8 + ...
                    randn*1e-10*x^7 + ...
                    randn*1e-10*x^6 + ...
                    randn*1e-10*x^5 +...
                    randn*1e-10*x^4 + ...
                    randn*1e-10*x^3 + ...
                    randn*1e-10*x^2 + ...
                    randn*1e-10*x^1 + ...
                    randn*1e-10; 

Z_perturbed = solve(WP_perturbed) %Hitta rötterna exakt till det störda polynomet 
plot(real(Z_perturbed),imag(Z_perturbed),'xr','linewidth',2) %Rita upp de störda rötterna
legend('Exact roots','Perturbed roots','location','northwest')

pause(0.1)

end

