function root = plot_poly(R,lmax,lmin)
B = zeros(3,3);
BB = R*R';
B(1,1) = R(1,1)^2 + R(1,2)^2 + R(1,3)^2;
B(2,2) = R(2,1)^2 + R(2,2)^2 + R(2,3)^2;
B(3,3) = R(3,1)^2 + R(3,2)^2 + R(3,3)^2;
B(1,2) = R(1,1)*R(2,1) + R(1,2)*R(2,2) + R(1,3)*R(2,3);
B(1,3) = R(1,1)*R(3,1) + R(1,2)*R(3,2) + R(1,3)*R(3,3);
B(2,3) = R(2,1)*R(3,1) + R(2,2)*R(3,2) + R(2,3)*R(3,3);
B(3,2) = B(2,3); B(3,1) = B(1,3); B(2,1) = B(1,2);
detR   = R(1,1)*(R(2,2)*R(3,3) - R(3,2)*R(2,3))...
       + R(1,2)*(R(2,3)*R(3,1) - R(2,1)*R(3,3))...
       + R(1,3)*(R(2,1)*R(3,2) - R(2,2)*R(3,1))
   %det_R = det(R)  % gives same answer
%lambda = eig(B);
%B-BB   % shows that B is indeed R^2, computing symmetric part once
%s = sqrt(lambda)*sign(detR);
%lam1 =  s(1)+s(2)+s(3)   % write out evalues based on SVD
%lam2 =  s(1)-s(2)-s(3)
%lam3 = -s(1)+s(2)-s(3)
%lam4 = -s(1)-s(2)+s(3)   % same as based on quaternions
p3 = -(B(1,1) + B(2,2) + B(3,3));
r = p3^2 + 4*(B(1,2)^2 + B(2,3)^2 + B(1,3)^2 -...
              B(2,2)*B(3,3)-B(3,3)*B(1,1)-B(1,1)*B(2,2));
q4 = -8*detR;
p6 = 2*p3;
x = linspace(lmin,lmax,101);
y = r + x.*(q4 + x.*(p6 + x.^2));
z =-q4^2/64 + x.*((p3^2-r)/4 + x.*(p3 + x));
%[z4,z3] = check_roots(p3/3,q4/4,r) % cubic/quartic agree
[root,n,delta] = muller(p6,q4,r,lmax,20,.00000000001)
z0 = zeros(101,1);
figure
plot(x,y,x,z0,x,z)