Roots of a cubic polynomial by solving a system with sums of powers of roots

%By Chukwunomso Agunwamba
%finds the roots of cubic polynomials by employing the sums of powers
%of their roots to construct and solve a system of equations.
%This is for cubics with distinct roots. When roots are repeated, a GCF is run
%using the derivative in order to remove the repeated root.

##This mathematical code is provided as is, for the end user to be able to use it
##and to edit it for educational purposes. However, any modified version must
##include my name as the original producer of the code.

##The name of the mfile code is findroots3. It has two input arguments.
##The first input is the polynomial coefficients vector in increasing powers.  
##The second input is an index to choose from one of the two possible linearly
##independent root interpolating polynomials the system solution calculates.
##This is because to find the coefficients of the interpolating polynomial,
##a quadratic equation was solved as a subsystem of the overall system of
##equations.
##
##There are two possible outputs. The first output contains the roots. The
##second output is the roots interpolating polynomial chosen by tq.
##The index, tq, is 1 or 2.
##
##
##Usage syntax example:
##
##
##findroots3(coeff([1 2 3]),1)
##ans =
##
##   2.000000000000025e+00 - 3.061439990403869e-13i
##   3.000000000000066e+00 - 5.334621633323877e-14i
##   9.999999999999438e-01 - 6.283862319378386e-14i

##[R, Phat] = findroots3(coeff([1 .2 .3]),1)
##R =
##
##   3.000000000000009e-01 - 2.473715676742927e-14i
##   9.999999999999842e-01 - 1.294103713078698e-14i
##   1.999999999999822e-01 + 1.044303532538038e-14i
##
##Phat =
##
## Columns 1 and 2:
##
##  -1.577192982456158e-01 + 1.880946403307244e-01i   6.921052631578980e-01 - 1.028595084845726e+00i
##
## Column 3:
##
##  -4.999999999999998e-01 + 8.660254037844392e-01i