function R = zroots(a,epsilon,N,p)
%
% R = zroots(a,epsilon,N,p)
%
% This is a successive application of Laguerre's Method applied to the
% polynomial with coefficients (with descending degree) in a and starting
% the iteration process at p.  If no p is specified, p=0 is assumed.  The
% process of finding each zero continues until either N iterations are
% performed or |q(p)|<epsilon.  An array of the (complex) roots is returned.
%

% Troy Henderson
% Assistant Professor
% University of Mobile
% Updated 7 February, 2008

ad = a;
m = length(a)-1;
for j=1:m-1
	r = laguerre(ad,epsilon,N);
	rrabs = abs(real(r));
	riabs = abs(imag(r));
	if (riabs < 2*epsilon*rrabs), r = real(r); end
	if (rrabs < 2*epsilon*riabs), r = imag(r); end
	R(j,1) = r;
	ad = deflate(ad,r);
end
R(m,1) = -ad(2)/ad(1);
for j=1:m
	R(j,1)=laguerre(a,epsilon,N,R(j,1));
end
R=sort(R);
endfunction
