disp(' cubic_poly_roots.m  ver 1.0  December 1, 2011');
disp(' by Tom Irvine');
disp(' ');
disp(' This scripts finds the roots of ');
disp(' ');
disp(' f(x)=ax^3+bx^2+cx+d ');
disp(' ');
%
clear x;
clear z;
clear w;
clear a;
clear b;
clear c;
clear d;
%
a=input(' enter a ');
b=input(' enter b ');
c=input(' enter c ');
d=input(' enter d ');
%
j=sqrt(-1);
%
p=(1/a)*((-b^2/(3*a))+c);
q=(1/(27*a^3))*((2*b^3)-(9*a*b*c)+(27*a^2*d));
%
r1=sqrt(q^2+4*p^3/27);
%
N1=(-q+r1)/2;
N2=(-q-r1)/2;
%
aa=real(N1);
bb=imag(N1);
%
alpha=(1/3)*atan2(bb,aa);
%
z(1)=((aa^2+bb^2)^(1/6))*(cos(alpha)+j*sin(alpha));
z(2)=z(1)*(-1+j*sqrt(3))/2;
z(3)=z(1)*(-1-j*sqrt(3))/2;
%
aa=real(N2);
bb=imag(N2);
%
alpha=(1/3)*atan2(bb,aa);
%
z(4)=((aa^2+bb^2)^(1/6))*(cos(alpha)+j*sin(alpha));
z(5)=z(4)*(-1+j*sqrt(3))/2;
z(6)=z(4)*(-1-j*sqrt(3))/2;
%
disp(' ');
disp(' roots ');
%
for(i=1:3)
    w(i)=z(i)-p/(3*z(i));
    x(i)=w(i)-b/(3*a);
    out1=sprintf(' %8.4g +j( %8.4g)',real(x(i)),imag(x(i)));
    disp(out1);
end
%
figure(1);
f=zeros(100,1);
L=max(abs(x));
U=min(abs(x));
P=max([L U]);
P=P*2;
r=linspace(-P,P,100);
for(i=1:100)
    f(i)=a*r(i)^3+b*r(i)^2+c*r(i)+d;
end
plot(r,f);
out1=sprintf(' f(x)= %dx^3 +%dx^2 +%dx +%d ',a,b,c,d);
title(out1);
xlabel('x');
grid on;
%