% JPD suggested this case
clear all;close all;clc;
c=[':.b'; ':sr'; ':og'; ':pm'; ':.k'];
A=[1 0 0 0; 0 2 0 0; 0 0 3 0; 0 0 0 4]; %place roots on diag
P=poly(A);
P1=polyder(P);
R=roots(P);
R1=roots(P1);
sprintf('Roots: %d, %d, %d, %d\n\n',R);
sprintf('Roots: %d, %d, %d, %d\n\n',R1);

y=[]; x=1+eps*j;
  while x~=0;
%clc;

%Now apply Newton method;
%Define starting value
%x0=input('starting value=');
%n=1;x(1)=x0+1e-6j;
x=6*(rand(1))  +3*(rand(1)-.5)*j;
%x=conj(x);

y=[];
for n=1:40
y(n)=x; %save sequence
x = x - 0.2*polyval(P,x)./polyval(P1,x);
%Plot results
%plot(y,c(rn,:),'MarkerSize',2); hold on; grid on;
%plot(1:length(R),R,'x','MarkerSize',5);hold on

  end %for

  rn=round(real(y(end)));  
  plot(y,c(rn,:),'MarkerSize',2); hold on; grid on;
  xlabel('\Re x'); ylabel('\Im x');
  xlim([0 5]);ylim([-2 2]);
  
%drawnow;%hold off;
disp(sprintf('y(end)=%g\n',real(y(end))));drawnow;
pause(.01); %input('Any key to continue');
  end %while
