Many physical phenomena can be described using partial differential equations (PDEs), they usually are non-linear and, due to the non-linearity, multiple different solutions can exist. Understanding this property and the relation between the different solutions and some parameters is crucial to predict the evolution of a system when the parameters can vary. Even if these phenomena can be represented with bifurcation diagrams, obtaining them analytically or theoretically is impossible for almost every interesting problem. For this reason, one would like to compute them numerically, but, due to the complexity of the task, several different techniques must be used together to perform it in an accurate and efficient way.