Solving quartic equation? (Cardano/Ferrari)

Question

today I've written a little Code-Snippet that is based upon the steps that are mentionned in this wikipedia-Article to solve a general quartic polynom.

Here's my matlab-implementation:

function ply4_roots = roots_cardano (dbl_p4, dbl_p3, dbl_p2, ...
    dbl_p1, dbl_p0)

    % solving a quartic equation with method of Ferrari

    dbl_alpha = -(3/8)*(dbl_p3^2)/(dbl_p4^2)+(dbl_p2/dbl_p4);
    dbl_beta = (1/8)*(dbl_p3^3)/(dbl_p4^3)-(1/2) ...
        *(dbl_p3*dbl_p2)/(dbl_p4^2)+ (dbl_p1/dbl_p4);
    dbl_gamma = -(3/256)*(dbl_p3^4)/(dbl_p4^4)+(1/16) ...
        *(dbl_p3^2*dbl_p2)/(dbl_p4^3) - (1/4) ...
        *(dbl_p3*dbl_p1)/(dbl_p4^2)+(dbl_p0/dbl_p4);

    dbl_pp = -(1/12)*dbl_alpha^2-dbl_gamma;
    dbl_qq = -(1/108)*dbl_alpha^3+(1/3)*dbl_alpha* ...
        dbl_gamma-(1/8)*dbl_beta^2;

    if dbl_pp~=0

        dbl_uu = (-(1/2)*dbl_qq+((1/4)*dbl_qq^2+(1/27)* ...
            dbl_pp^3)^(1/2))^(1/3);
        dbl_yy = -(5/6)*dbl_alpha+dbl_uu-(1/3)*(dbl_pp/dbl_uu);

    else

        dbl_yy = -(5/6)*dbl_alpha-(dbl_qq)^(1/3);

    end

    dbl_ww = (dbl_alpha+2*dbl_yy)^(1/2);

    ply4_roots(1) = -(1/4)*(dbl_p3/dbl_p4)+(1/2)*(1*dbl_ww+ ...
        1*(-(dbl_alpha+2*dbl_yy)-2*(dbl_alpha+ ...
        1*(dbl_beta/dbl_ww)))^(1/2));

    ply4_roots(2) = -(1/4)*(dbl_p3/dbl_p4)+(1/2)*(1*dbl_ww+ ...
        (-1)*(-(dbl_alpha+2*dbl_yy)-2*(dbl_alpha+ ...
        1*(dbl_beta/dbl_ww)))^(1/2));

    ply4_roots(3) = -(1/4)*(dbl_p3/dbl_p4)+(1/2)*((-1)*dbl_ww+ ...
        (-1)*(-(dbl_alpha+2*dbl_yy)-2*(dbl_alpha+ ...
        (-1)*(dbl_beta/dbl_ww)))^(1/2));

    ply4_roots(4) = -(1/4)*(dbl_p3/dbl_p4)+(1/2)*((-1)*dbl_ww+ ...
        1*(-(dbl_alpha+2*dbl_yy)-2*(dbl_alpha+ ...
        (-1)*(dbl_beta/dbl_ww)))^(1/2));

    % sort results with bubble-sort
    for ii=4:-1:1

        for idx=1:3

            if ply4_roots(idx)>ply4_roots(idx+1)

                value=ply4_roots(idx);
                ply4_roots(idx) = ply4_roots(idx+1);
                ply4_roots(idx+1) = value;

            end

        end

    end

end

In principle, it seems that the script works for examples like:

>> roots_cardano(1,-4,-1,16,-12)

ans =

  -2     1     2     3

But, if I try it for

>> roots_cardano(2, -4, -3, 7, -2)

ans =

  -1.3660 + 0.0000i   0.3660 - 0.0000i   1.0000 + 0.0000i   2.0000 - 0.0000i

I get "some" complex results which shouldn't be there. Matlab's root-finder gives:

>> roots([2 -4 -3 7 -2])

ans =

   -1.3660
    2.0000
    1.0000
    0.3660

what underlines that there are no complex parts in the solution.

For another example, it is somehow 'reversed':

>> roots_cardano(1,-14,44,-40,0)

ans =

         0    2.0000    2.0000   10.0000

... and matlab gives:

>> roots([1 -14 44 -40 0])

ans =

        0          
  10.0000          
   2.0000 + 0.0000i
   2.0000 - 0.0000i

and one last example:

    coef_cply4 = 2.03975258797202e-14
    coef_cply3 = 5.29694855354049e-14
    coef_cply2 = 2.83428012361253e-14
    coef_cply1 = 3.02190058771691e-15
    coef_cply0 = 2.80521515614180e-20

    >> roots_cardano(coef_cply4, coef_cply3, coef_cply2, coef_cply1, coef_cply0)

    ans =

      -1.9100 + 0.0000i  -0.5444 - 0.0000i  -0.1425 + 0.0000i  -0.0000 - 0.0000i

>> roots([coef_cply4 coef_cply3 coef_cply2 coef_cply1 coef_cply0])

ans =

   -1.9100
   -0.5444
   -0.1425
   -0.0000

... so I guess that there could be somewhere a mistake in my script. Unfortunately I did not find it until now.

It would be great if someone could help me ...

Thank you in advance!