function [r1,r2]=findroots(coeff,r0)
% [r1,r2]=findroots(coeff,r0)
% determines numerically, using Newton iterations, 
% the two zeros r1 and r2 of the radii polynomial 
% coeff(1)*r^N + ... + coeff(N+1)*r^0, where N=length(coeff)-1
% using the initial guesses r0=[r0(1) r0(2)]

maxsteps=40;
tol=1e-14;

nsteps=0;
solfound=0;
del=tol+1;
r=r0(1);

n=length(coeff);
pol=@(r) sum(coeff.*r.^(n-1:-1:0));
dpol=@(r) sum(coeff(1:n-1).*r.^(n-2:-1:0).*(n-1:-1:1));

while solfound==0 && nsteps<maxsteps;
    if del>tol
      dr=-pol(r)/dpol(r);
      del=norm(dr);
      r=r+dr;
      nsteps=nsteps+1;
    else
%        disp('Newton converged')
        solfound=1;
        r1=r;
    end    
end

if solfound
    nsteps=0;
    solfound=0;
    del=tol+1;
    r=r0(2);
    while solfound==0 && nsteps<maxsteps;
      if del>tol
        dr=-pol(r)/dpol(r);
        del=norm(dr);
        r=r+dr;
        nsteps=nsteps+1;
      else
%        disp('Newton converged')
        solfound=1;
        r2=r;
      end    
    end
end

if ~solfound
    r1=NaN;
    r2=NaN;
end

end
