The Muller's method is supposed for finding the real solution. Your code is attempting to find the complex solution. I am wondering if it has any other side-effect.
The Muller's method is based on secant method (or in a similar way) to iteratively find the solution. It cannot guarantee to find all solutions of a functions. It just find a solution (normally real) within the specified interval.
For you code, you may find another solution by change the sign of the following line
D = -(b.^2 - 4*f(p2).*d).^(1/2);