What is a simple way to find real roots of a (cubic) polynomial?

For a cubic polynomial there are closed form solutions, but they are not particularly well suited for numerical calculus.

I'd do the following for the cubic case: any cubic polynomial has at least one real root, you can find it easily with Newton's method. Then, you use deflation to get the remaining quadratic polynomial to solve, see my answer there for how to do this latter step correctly.

One word of caution: if the discriminant is close to zero, there will be a numerically multiple real root, and Newton's method will miserably fail. Moreover, since to the vicinity of the root, the polynomial is like (x - x0)^2, you'll lose half your significant digits (since P(x) will be < epsilon as soon as x - x0 < sqrt(epsilon)). So you may want to rule this out and use the closed form solution in this particular case, or solve the derivative polynomial.

If you want to find roots in a given interval, check Sturm's theorem.

A more general (complex) algorithm for generic polynomial solving is Jenkins-Traub algorithm. This is clearly overkill here, but it works well on cubics. Usually, you use a third-party implementation.

Since you do C, using the GSL is surely your best bet.

Another generic method is to find the eigenvalues of the companion matrix with eg. balanced QR decomposition, or reduction to Householder form. This is the approach taken by GSL.

If you don't want to use the closed from solutions (or expect polynoms of larger order), the most obvious method would be to calculate approximate roots by using Newton's method.

Unfortunately it's not possible to decide which roots you will get when iterating, although it depends on the starting value.

Also see here.

See Solving quartics and cubics for graphics by D Herbison-Evans, published in Graphics Gems V.