R code to solve Cubic Equation

This is create function to solve cubic equation.

note: default form of a cubic equation in this "cubicEq" function is: a*X^3+b*X^2+c*X+d=0. 

cubicEq<-function(a,b,c,d)

{  

   x1<-c()

   x2<-c()

   x3<-c()

   for ( i in 1:2)

   {

      f<-(3*c[i]/a[i]-(b[i]^2/(a[i]^2)))/3

      g<-(2*b[i]^3/a[i]^3-(9*b[i]*c[i]/a[i]^2)+27*d[i]/a[i])/27

      h<-(g^2/4)+(f^3/27)

      if (h<=0)

         {

             count<-(g*g/4-h)^(1/2)

             j<-count^(1/3)

             k<-acos(-g/(2*count))

             L<-j*(-1)

             M<-cos(k/3)

             N<-3^(1/2)*sin(k/3)

             p<-(-1)*(b[i]/(3*a[i]))

             x1[i]<-2*j*cos(k/3)-b[i]/(3*a[i])

             x2[i]<-L*(M+N)+p

             x3[i]<-L*(M-N)+p

          }

       if ( h>0 )

          {

             R<- -g/2+h^(1/2)

             S<-R^(1/3)

             T<- -g/2-h^(1/2)

             U<--(-T)^(1/3)

             x1[i]<-(S+U)-(b[i]/(3*a[i]))

             x2[i]<- complex(real=-(S+U)/2-b/(3*a[i]),imaginary=((S-U)*sqrt(3)/2))

             x3[i]<- complex(real=-(S+U)/2-b/(3*a[i]),imaginary=-((S-U)*sqrt(3)/2))

          }

      }

return(list(x1=x1,x2=x2,x3=x3))

}