What am I trying to do? I'm trying to implement a simulation of a chaotic billiard system, following the algorithm in this excerpt. How am I trying it? Using numpy and matplotlib, I implemented the following code What is the problem? When calculating phi_new, the equation has two solutions (assuming the boundary is convex, which is.) I must enforce that phi_new is the solution which is different from phi but I don't know how to do that. Are there more issues with the code? What should the output be? A phase space diagram of S x Alpha, looking like this. Any help is very appreciated! Thanks in advance.