I am trying to implement a root finding algorithm. I am using the hybrid Newton-Raphson algorithm found in numerical recipes that works pretty nicely. But I have a problem in bracketing the root.

While implementing the root finding algorithm I realised that in several cases my functions have 1 real root and all the other imaginary (several of them, usually 6 or 9). The only root I am interested is in the real one so the problem is not there. The thing is that the function approaches the root like a cubic function, touching with the point the y=0 axis...

Newton-Rapson method needs some brackets of different sign and all the bracketing methods I found don't work for this specific case.

What can I do? It is pretty important to find that root in my program...

EDIT: more problems: sometimes due to reaaaaaally small numerical errors, say a variation of 1e-6 in some value the "cubic" function does NOT have that real root, it is just imaginary with a neglectable imaginary part... (checked with matlab)

EDIT 2: Much more information about the problem.

Ok, I need root finding algorithm.

Info I have:

• The root I need to find is between [0-1] , if there are more roots outside that part I am not interested in them.
• The root is real, there may be imaginary roots, but I don't want them.
• Probably all the rest of the roots will be imaginary
• The root may be double in that point, but I think that actually doesn't mater in numerical analysis problems
• I need to use the root finding algorithm several times during the overall calculations, but the function will always be a polynomial
• In one of the particular cases of the root finding, my polynomial will be similar to a quadratic function that touches Y=0 with the point. Example of a real case:
• The coefficient may not be 100% precise and that really slight imprecision may make the function not to touch the Y=0 axis.
• I cannot solve for this specific case because in other cases it may be that the polynomial is pretty normal and doesn't make any "strange" thing.
• The method I am actually using is NewtonRaphson hybrid, where if the derivative is really small it makes a bisection instead of NewRaph (found in numerical recipes).

Matlab's answer to the function on the image: roots:

0.853553390593276 + 0.353553390593278i
0.853553390593276 - 0.353553390593278i
0.146446609406726 + 0.353553390593273i
0.146446609406726 - 0.353553390593273i
0.499999999999996 + 0.000000040142134i
0.499999999999996 - 0.000000040142134i

The function is a real example I prepared where I know that the answer I want is 0.5

Note: I still haven't check completely some of the answers I you people have give me (Thank you!), I am just trying to give al the information I already have to complete the question.