function z = poly_root(P)

% poly_root - find the roots of a polynomial
%
%   z = root_low_poly(C);
%
%   computes the roots of the polynomial whose coefficients
%   are the elements of the vector C. If C has N+1 components,
%   the polynomial is C(1)*X^N + ... + C(N)*X + C(N+1).
%
%   Same input structure as roots().
%
%   Special code for degrees <= 4.
%   This method can take a list of polynomials, 
%   i.e. C(:,i) is the coefficient of the ith polynomial
%   (fast parallel resolution).
%
%
%   Copyright (c) 2011 Gabriel Peyré

if size(P,1)==1 && size(P,2)~=1
    P = P';
end

if size(P,1)==2
    P = P(end:-1:1,:);    
    z = -P(1,:)./P(2,:);
    
elseif size(P,1)==3
    P = P(end:-1:1,:);
    
    delta = P(2,:).^2 - 4*P(1,:).*P(3,:);
    z = [(-P(2,:)+sqrt(delta))./( 2*P(3,:) ); (-P(2,:)-sqrt(delta))./( 2*P(3,:) ) ];
    
elseif size(P,1)==4    
    P = P(end:-1:1,:);
   
   % normalize polynomial
    for i=1:3
        P(i,:) = P(i,:)./P(4,:);
    end
    P(4,:) = 1;
    
    p = -P(3,:).*P(3,:)/9+P(2,:)/3;
    q = 2*P(3,:).^3/27-P(3,:).*P(2,:)/3+P(1,:);
    S = [-p.^3; q; P(4,:)];
    rr = poly_root(S(end:-1:1,:));   % remember P(4,:) = 1
    u = rr(1,:).^(1/3);
    v = rr(2,:).^(1/3);
    z = [u+v-P(3,:)/3;
        -0.5*(u+v)-P(3,:)/3 + 0.5i*sqrt(3).*(u-v);
        -0.5*(u+v)-P(3,:)/3 - 0.5i*sqrt(3).*(u-v)];
    
    
elseif size(P,1)==5
    
    a = P(1,:);
    b = P(2,:);
    c = P(3,:);
    d = P(4,:);
    e = P(5,:);
    z(1,:) = -(1./4).*b./a+(1./2).*sqrt((1./4).*b.^2./a.^2-(2./3).*c./a+(1./6).*(-36.*d.*c.*b-288.*e.*c.*a+108.*d.^2.*a+108.*e.*b.^2+8.*c.^3+12.*sqrt(-3.*d.^2.*b.^2.*c.^2+12.*e.*b.^2.*c.^3-54.*d.*c.*b.^3.*e+12.*d.^3.*b.^3+81.*e.^2.*b.^4+12.*d.^2.*a.*c.^3+384.*a.^2.*e.^2.*c.^2-48.*a.*e.*c.^4+240.*d.*b.*a.*e.*c.^2+18.*d.^2.*b.^2.*a.*e+576.*d.*b.*e.^2.*a.^2-54.*d.^3.*c.*b.*a-432.*e.*c.*a.^2.*d.^2-432.*e.^2.*c.*a.*b.^2-768.*a.^3.*e.^3+81.*d.^4.*a.^2)).^(1./3)./a+(2./3).*(-3.*d.*b+12.*a.*e+c.^2)./(a.*(-36.*d.*c.*b-288.*e.*c.*a+108.*d.^2.*a+108.*e.*b.^2+8.*c.^3+12.*sqrt(-3.*d.^2.*b.^2.*c.^2+12.*e.*b.^2.*c.^3-54.*d.*c.*b.^3.*e+12.*d.^3.*b.^3+81.*e.^2.*b.^4+12.*d.^2.*a.*c.^3+384.*a.^2.*e.^2.*c.^2-48.*a.*e.*c.^4+240.*d.*b.*a.*e.*c.^2+18.*d.^2.*b.^2.*a.*e+576.*d.*b.*e.^2.*a.^2-54.*d.^3.*c.*b.*a-432.*e.*c.*a.^2.*d.^2-432.*e.^2.*c.*a.*b.^2-768.*a.^3.*e.^3+81.*d.^4.*a.^2)).^(1./3)))+(1./2).*sqrt((1./2).*b.^2./a.^2-(4./3).*c./a-(1./6).*(-36.*d.*c.*b-288.*e.*c.*a+108.*d.^2.*a+108.*e.*b.^2+8.*c.^3+12.*sqrt(-3.*d.^2.*b.^2.*c.^2+12.*e.*b.^2.*c.^3-54.*d.*c.*b.^3.*e+12.*d.^3.*b.^3+81.*e.^2.*b.^4+12.*d.^2.*a.*c.^3+384.*a.^2.*e.^2.*c.^2-48.*a.*e.*c.^4+240.*d.*b.*a.*e.*c.^2+18.*d.^2.*b.^2.*a.*e+576.*d.*b.*e.^2.*a.^2-54.*d.^3.*c.*b.*a-432.*e.*c.*a.^2.*d.^2-432.*e.^2.*c.*a.*b.^2-768.*a.^3.*e.^3+81.*d.^4.*a.^2)).^(1./3)./a-(2./3).*(-3.*d.*b+12.*a.*e+c.^2)./(a.*(-36.*d.*c.*b-288.*e.*c.*a+108.*d.^2.*a+108.*e.*b.^2+8.*c.^3+12.*sqrt(-3.*d.^2.*b.^2.*c.^2+12.*e.*b.^2.*c.^3-54.*d.*c.*b.^3.*e+12.*d.^3.*b.^3+81.*e.^2.*b.^4+12.*d.^2.*a.*c.^3+384.*a.^2.*e.^2.*c.^2-48.*a.*e.*c.^4+240.*d.*b.*a.*e.*c.^2+18.*d.^2.*b.^2.*a.*e+576.*d.*b.*e.^2.*a.^2-54.*d.^3.*c.*b.*a-432.*e.*c.*a.^2.*d.^2-432.*e.^2.*c.*a.*b.^2-768.*a.^3.*e.^3+81.*d.^4.*a.^2)).^(1./3))+(c.*b./a.^2-2.*d./a-(1./4).*b.^3./a.^3)./sqrt((1./4).*b.^2./a.^2-(2./3).*c./a+(1./6).*(-36.*d.*c.*b-288.*e.*c.*a+108.*d.^2.*a+108.*e.*b.^2+8.*c.^3+12.*sqrt(-3.*d.^2.*b.^2.*c.^2+12.*e.*b.^2.*c.^3-54.*d.*c.*b.^3.*e+12.*d.^3.*b.^3+81.*e.^2.*b.^4+12.*d.^2.*a.*c.^3+384.*a.^2.*e.^2.*c.^2-48.*a.*e.*c.^4+240.*d.*b.*a.*e.*c.^2+18.*d.^2.*b.^2.*a.*e+576.*d.*b.*e.^2.*a.^2-54.*d.^3.*c.*b.*a-432.*e.*c.*a.^2.*d.^2-432.*e.^2.*c.*a.*b.^2-768.*a.^3.*e.^3+81.*d.^4.*a.^2)).^(1./3)./a+(2./3).*(-3.*d.*b+12.*a.*e+c.^2)./(a.*(-36.*d.*c.*b-288.*e.*c.*a+108.*d.^2.*a+108.*e.*b.^2+8.*c.^3+12.*sqrt(-3.*d.^2.*b.^2.*c.^2+12.*e.*b.^2.*c.^3-54.*d.*c.*b.^3.*e+12.*d.^3.*b.^3+81.*e.^2.*b.^4+12.*d.^2.*a.*c.^3+384.*a.^2.*e.^2.*c.^2-48.*a.*e.*c.^4+240.*d.*b.*a.*e.*c.^2+18.*d.^2.*b.^2.*a.*e+576.*d.*b.*e.^2.*a.^2-54.*d.^3.*c.*b.*a-432.*e.*c.*a.^2.*d.^2-432.*e.^2.*c.*a.*b.^2-768.*a.^3.*e.^3+81.*d.^4.*a.^2)).^(1./3)))); 
    z(2,:) = -(1./4).*b./a+(1./2).*sqrt((1./4).*b.^2./a.^2-(2./3).*c./a+(1./6).*(-36.*d.*c.*b-288.*e.*c.*a+108.*d.^2.*a+108.*e.*b.^2+8.*c.^3+12.*sqrt(-3.*d.^2.*b.^2.*c.^2+12.*e.*b.^2.*c.^3-54.*d.*c.*b.^3.*e+12.*d.^3.*b.^3+81.*e.^2.*b.^4+12.*d.^2.*a.*c.^3+384.*a.^2.*e.^2.*c.^2-48.*a.*e.*c.^4+240.*d.*b.*a.*e.*c.^2+18.*d.^2.*b.^2.*a.*e+576.*d.*b.*e.^2.*a.^2-54.*d.^3.*c.*b.*a-432.*e.*c.*a.^2.*d.^2-432.*e.^2.*c.*a.*b.^2-768.*a.^3.*e.^3+81.*d.^4.*a.^2)).^(1./3)./a+(2./3).*(-3.*d.*b+12.*a.*e+c.^2)./(a.*(-36.*d.*c.*b-288.*e.*c.*a+108.*d.^2.*a+108.*e.*b.^2+8.*c.^3+12.*sqrt(-3.*d.^2.*b.^2.*c.^2+12.*e.*b.^2.*c.^3-54.*d.*c.*b.^3.*e+12.*d.^3.*b.^3+81.*e.^2.*b.^4+12.*d.^2.*a.*c.^3+384.*a.^2.*e.^2.*c.^2-48.*a.*e.*c.^4+240.*d.*b.*a.*e.*c.^2+18.*d.^2.*b.^2.*a.*e+576.*d.*b.*e.^2.*a.^2-54.*d.^3.*c.*b.*a-432.*e.*c.*a.^2.*d.^2-432.*e.^2.*c.*a.*b.^2-768.*a.^3.*e.^3+81.*d.^4.*a.^2)).^(1./3)))-(1./2).*sqrt((1./2).*b.^2./a.^2-(4./3).*c./a-(1./6).*(-36.*d.*c.*b-288.*e.*c.*a+108.*d.^2.*a+108.*e.*b.^2+8.*c.^3+12.*sqrt(-3.*d.^2.*b.^2.*c.^2+12.*e.*b.^2.*c.^3-54.*d.*c.*b.^3.*e+12.*d.^3.*b.^3+81.*e.^2.*b.^4+12.*d.^2.*a.*c.^3+384.*a.^2.*e.^2.*c.^2-48.*a.*e.*c.^4+240.*d.*b.*a.*e.*c.^2+18.*d.^2.*b.^2.*a.*e+576.*d.*b.*e.^2.*a.^2-54.*d.^3.*c.*b.*a-432.*e.*c.*a.^2.*d.^2-432.*e.^2.*c.*a.*b.^2-768.*a.^3.*e.^3+81.*d.^4.*a.^2)).^(1./3)./a-(2./3).*(-3.*d.*b+12.*a.*e+c.^2)./(a.*(-36.*d.*c.*b-288.*e.*c.*a+108.*d.^2.*a+108.*e.*b.^2+8.*c.^3+12.*sqrt(-3.*d.^2.*b.^2.*c.^2+12.*e.*b.^2.*c.^3-54.*d.*c.*b.^3.*e+12.*d.^3.*b.^3+81.*e.^2.*b.^4+12.*d.^2.*a.*c.^3+384.*a.^2.*e.^2.*c.^2-48.*a.*e.*c.^4+240.*d.*b.*a.*e.*c.^2+18.*d.^2.*b.^2.*a.*e+576.*d.*b.*e.^2.*a.^2-54.*d.^3.*c.*b.*a-432.*e.*c.*a.^2.*d.^2-432.*e.^2.*c.*a.*b.^2-768.*a.^3.*e.^3+81.*d.^4.*a.^2)).^(1./3))+(c.*b./a.^2-2.*d./a-(1./4).*b.^3./a.^3)./sqrt((1./4).*b.^2./a.^2-(2./3).*c./a+(1./6).*(-36.*d.*c.*b-288.*e.*c.*a+108.*d.^2.*a+108.*e.*b.^2+8.*c.^3+12.*sqrt(-3.*d.^2.*b.^2.*c.^2+12.*e.*b.^2.*c.^3-54.*d.*c.*b.^3.*e+12.*d.^3.*b.^3+81.*e.^2.*b.^4+12.*d.^2.*a.*c.^3+384.*a.^2.*e.^2.*c.^2-48.*a.*e.*c.^4+240.*d.*b.*a.*e.*c.^2+18.*d.^2.*b.^2.*a.*e+576.*d.*b.*e.^2.*a.^2-54.*d.^3.*c.*b.*a-432.*e.*c.*a.^2.*d.^2-432.*e.^2.*c.*a.*b.^2-768.*a.^3.*e.^3+81.*d.^4.*a.^2)).^(1./3)./a+(2./3).*(-3.*d.*b+12.*a.*e+c.^2)./(a.*(-36.*d.*c.*b-288.*e.*c.*a+108.*d.^2.*a+108.*e.*b.^2+8.*c.^3+12.*sqrt(-3.*d.^2.*b.^2.*c.^2+12.*e.*b.^2.*c.^3-54.*d.*c.*b.^3.*e+12.*d.^3.*b.^3+81.*e.^2.*b.^4+12.*d.^2.*a.*c.^3+384.*a.^2.*e.^2.*c.^2-48.*a.*e.*c.^4+240.*d.*b.*a.*e.*c.^2+18.*d.^2.*b.^2.*a.*e+576.*d.*b.*e.^2.*a.^2-54.*d.^3.*c.*b.*a-432.*e.*c.*a.^2.*d.^2-432.*e.^2.*c.*a.*b.^2-768.*a.^3.*e.^3+81.*d.^4.*a.^2)).^(1./3)))); 
    z(3,:) = -(1./4).*b./a-(1./2).*sqrt((1./4).*b.^2./a.^2-(2./3).*c./a+(1./6).*(-36.*d.*c.*b-288.*e.*c.*a+108.*d.^2.*a+108.*e.*b.^2+8.*c.^3+12.*sqrt(-3.*d.^2.*b.^2.*c.^2+12.*e.*b.^2.*c.^3-54.*d.*c.*b.^3.*e+12.*d.^3.*b.^3+81.*e.^2.*b.^4+12.*d.^2.*a.*c.^3+384.*a.^2.*e.^2.*c.^2-48.*a.*e.*c.^4+240.*d.*b.*a.*e.*c.^2+18.*d.^2.*b.^2.*a.*e+576.*d.*b.*e.^2.*a.^2-54.*d.^3.*c.*b.*a-432.*e.*c.*a.^2.*d.^2-432.*e.^2.*c.*a.*b.^2-768.*a.^3.*e.^3+81.*d.^4.*a.^2)).^(1./3)./a+(2./3).*(-3.*d.*b+12.*a.*e+c.^2)./(a.*(-36.*d.*c.*b-288.*e.*c.*a+108.*d.^2.*a+108.*e.*b.^2+8.*c.^3+12.*sqrt(-3.*d.^2.*b.^2.*c.^2+12.*e.*b.^2.*c.^3-54.*d.*c.*b.^3.*e+12.*d.^3.*b.^3+81.*e.^2.*b.^4+12.*d.^2.*a.*c.^3+384.*a.^2.*e.^2.*c.^2-48.*a.*e.*c.^4+240.*d.*b.*a.*e.*c.^2+18.*d.^2.*b.^2.*a.*e+576.*d.*b.*e.^2.*a.^2-54.*d.^3.*c.*b.*a-432.*e.*c.*a.^2.*d.^2-432.*e.^2.*c.*a.*b.^2-768.*a.^3.*e.^3+81.*d.^4.*a.^2)).^(1./3)))+(1./2).*sqrt((1./2).*b.^2./a.^2-(4./3).*c./a-(1./6).*(-36.*d.*c.*b-288.*e.*c.*a+108.*d.^2.*a+108.*e.*b.^2+8.*c.^3+12.*sqrt(-3.*d.^2.*b.^2.*c.^2+12.*e.*b.^2.*c.^3-54.*d.*c.*b.^3.*e+12.*d.^3.*b.^3+81.*e.^2.*b.^4+12.*d.^2.*a.*c.^3+384.*a.^2.*e.^2.*c.^2-48.*a.*e.*c.^4+240.*d.*b.*a.*e.*c.^2+18.*d.^2.*b.^2.*a.*e+576.*d.*b.*e.^2.*a.^2-54.*d.^3.*c.*b.*a-432.*e.*c.*a.^2.*d.^2-432.*e.^2.*c.*a.*b.^2-768.*a.^3.*e.^3+81.*d.^4.*a.^2)).^(1./3)./a-(2./3).*(-3.*d.*b+12.*a.*e+c.^2)./(a.*(-36.*d.*c.*b-288.*e.*c.*a+108.*d.^2.*a+108.*e.*b.^2+8.*c.^3+12.*sqrt(-3.*d.^2.*b.^2.*c.^2+12.*e.*b.^2.*c.^3-54.*d.*c.*b.^3.*e+12.*d.^3.*b.^3+81.*e.^2.*b.^4+12.*d.^2.*a.*c.^3+384.*a.^2.*e.^2.*c.^2-48.*a.*e.*c.^4+240.*d.*b.*a.*e.*c.^2+18.*d.^2.*b.^2.*a.*e+576.*d.*b.*e.^2.*a.^2-54.*d.^3.*c.*b.*a-432.*e.*c.*a.^2.*d.^2-432.*e.^2.*c.*a.*b.^2-768.*a.^3.*e.^3+81.*d.^4.*a.^2)).^(1./3))-(c.*b./a.^2-2.*d./a-(1./4).*b.^3./a.^3)./sqrt((1./4).*b.^2./a.^2-(2./3).*c./a+(1./6).*(-36.*d.*c.*b-288.*e.*c.*a+108.*d.^2.*a+108.*e.*b.^2+8.*c.^3+12.*sqrt(-3.*d.^2.*b.^2.*c.^2+12.*e.*b.^2.*c.^3-54.*d.*c.*b.^3.*e+12.*d.^3.*b.^3+81.*e.^2.*b.^4+12.*d.^2.*a.*c.^3+384.*a.^2.*e.^2.*c.^2-48.*a.*e.*c.^4+240.*d.*b.*a.*e.*c.^2+18.*d.^2.*b.^2.*a.*e+576.*d.*b.*e.^2.*a.^2-54.*d.^3.*c.*b.*a-432.*e.*c.*a.^2.*d.^2-432.*e.^2.*c.*a.*b.^2-768.*a.^3.*e.^3+81.*d.^4.*a.^2)).^(1./3)./a+(2./3).*(-3.*d.*b+12.*a.*e+c.^2)./(a.*(-36.*d.*c.*b-288.*e.*c.*a+108.*d.^2.*a+108.*e.*b.^2+8.*c.^3+12.*sqrt(-3.*d.^2.*b.^2.*c.^2+12.*e.*b.^2.*c.^3-54.*d.*c.*b.^3.*e+12.*d.^3.*b.^3+81.*e.^2.*b.^4+12.*d.^2.*a.*c.^3+384.*a.^2.*e.^2.*c.^2-48.*a.*e.*c.^4+240.*d.*b.*a.*e.*c.^2+18.*d.^2.*b.^2.*a.*e+576.*d.*b.*e.^2.*a.^2-54.*d.^3.*c.*b.*a-432.*e.*c.*a.^2.*d.^2-432.*e.^2.*c.*a.*b.^2-768.*a.^3.*e.^3+81.*d.^4.*a.^2)).^(1./3)))); 
    z(4,:) = -(1./4).*b./a-(1./2).*sqrt((1./4).*b.^2./a.^2-(2./3).*c./a+(1./6).*(-36.*d.*c.*b-288.*e.*c.*a+108.*d.^2.*a+108.*e.*b.^2+8.*c.^3+12.*sqrt(-3.*d.^2.*b.^2.*c.^2+12.*e.*b.^2.*c.^3-54.*d.*c.*b.^3.*e+12.*d.^3.*b.^3+81.*e.^2.*b.^4+12.*d.^2.*a.*c.^3+384.*a.^2.*e.^2.*c.^2-48.*a.*e.*c.^4+240.*d.*b.*a.*e.*c.^2+18.*d.^2.*b.^2.*a.*e+576.*d.*b.*e.^2.*a.^2-54.*d.^3.*c.*b.*a-432.*e.*c.*a.^2.*d.^2-432.*e.^2.*c.*a.*b.^2-768.*a.^3.*e.^3+81.*d.^4.*a.^2)).^(1./3)./a+(2./3).*(-3.*d.*b+12.*a.*e+c.^2)./(a.*(-36.*d.*c.*b-288.*e.*c.*a+108.*d.^2.*a+108.*e.*b.^2+8.*c.^3+12.*sqrt(-3.*d.^2.*b.^2.*c.^2+12.*e.*b.^2.*c.^3-54.*d.*c.*b.^3.*e+12.*d.^3.*b.^3+81.*e.^2.*b.^4+12.*d.^2.*a.*c.^3+384.*a.^2.*e.^2.*c.^2-48.*a.*e.*c.^4+240.*d.*b.*a.*e.*c.^2+18.*d.^2.*b.^2.*a.*e+576.*d.*b.*e.^2.*a.^2-54.*d.^3.*c.*b.*a-432.*e.*c.*a.^2.*d.^2-432.*e.^2.*c.*a.*b.^2-768.*a.^3.*e.^3+81.*d.^4.*a.^2)).^(1./3)))-(1./2).*sqrt((1./2).*b.^2./a.^2-(4./3).*c./a-(1./6).*(-36.*d.*c.*b-288.*e.*c.*a+108.*d.^2.*a+108.*e.*b.^2+8.*c.^3+12.*sqrt(-3.*d.^2.*b.^2.*c.^2+12.*e.*b.^2.*c.^3-54.*d.*c.*b.^3.*e+12.*d.^3.*b.^3+81.*e.^2.*b.^4+12.*d.^2.*a.*c.^3+384.*a.^2.*e.^2.*c.^2-48.*a.*e.*c.^4+240.*d.*b.*a.*e.*c.^2+18.*d.^2.*b.^2.*a.*e+576.*d.*b.*e.^2.*a.^2-54.*d.^3.*c.*b.*a-432.*e.*c.*a.^2.*d.^2-432.*e.^2.*c.*a.*b.^2-768.*a.^3.*e.^3+81.*d.^4.*a.^2)).^(1./3)./a-(2./3).*(-3.*d.*b+12.*a.*e+c.^2)./(a.*(-36.*d.*c.*b-288.*e.*c.*a+108.*d.^2.*a+108.*e.*b.^2+8.*c.^3+12.*sqrt(-3.*d.^2.*b.^2.*c.^2+12.*e.*b.^2.*c.^3-54.*d.*c.*b.^3.*e+12.*d.^3.*b.^3+81.*e.^2.*b.^4+12.*d.^2.*a.*c.^3+384.*a.^2.*e.^2.*c.^2-48.*a.*e.*c.^4+240.*d.*b.*a.*e.*c.^2+18.*d.^2.*b.^2.*a.*e+576.*d.*b.*e.^2.*a.^2-54.*d.^3.*c.*b.*a-432.*e.*c.*a.^2.*d.^2-432.*e.^2.*c.*a.*b.^2-768.*a.^3.*e.^3+81.*d.^4.*a.^2)).^(1./3))-(c.*b./a.^2-2.*d./a-(1./4).*b.^3./a.^3)./sqrt((1./4).*b.^2./a.^2-(2./3).*c./a+(1./6).*(-36.*d.*c.*b-288.*e.*c.*a+108.*d.^2.*a+108.*e.*b.^2+8.*c.^3+12.*sqrt(-3.*d.^2.*b.^2.*c.^2+12.*e.*b.^2.*c.^3-54.*d.*c.*b.^3.*e+12.*d.^3.*b.^3+81.*e.^2.*b.^4+12.*d.^2.*a.*c.^3+384.*a.^2.*e.^2.*c.^2-48.*a.*e.*c.^4+240.*d.*b.*a.*e.*c.^2+18.*d.^2.*b.^2.*a.*e+576.*d.*b.*e.^2.*a.^2-54.*d.^3.*c.*b.*a-432.*e.*c.*a.^2.*d.^2-432.*e.^2.*c.*a.*b.^2-768.*a.^3.*e.^3+81.*d.^4.*a.^2)).^(1./3)./a+(2./3).*(-3.*d.*b+12.*a.*e+c.^2)./(a.*(-36.*d.*c.*b-288.*e.*c.*a+108.*d.^2.*a+108.*e.*b.^2+8.*c.^3+12.*sqrt(-3.*d.^2.*b.^2.*c.^2+12.*e.*b.^2.*c.^3-54.*d.*c.*b.^3.*e+12.*d.^3.*b.^3+81.*e.^2.*b.^4+12.*d.^2.*a.*c.^3+384.*a.^2.*e.^2.*c.^2-48.*a.*e.*c.^4+240.*d.*b.*a.*e.*c.^2+18.*d.^2.*b.^2.*a.*e+576.*d.*b.*e.^2.*a.^2-54.*d.^3.*c.*b.*a-432.*e.*c.*a.^2.*d.^2-432.*e.^2.*c.*a.*b.^2-768.*a.^3.*e.^3+81.*d.^4.*a.^2)).^(1./3))));
    
    % round very small imaginary entries
    I = find(abs(imag(z))<1e-6);
    z(I) = real(z(I));

else
    
    warning('Works only for degree <=4, using roots');
    z = zeros(size(P,1)-1, size(P,2));
    for i=1:size(P,2)
        z(:,i) = roots(P(:,i));
    end
    
end


% I = find( imag(z)~=0 );
% z(I) = NaN;