# Optimal way to find stationary solutions of the PDE

I am researching heat diffusion in an optical element irradiated by laser. This problem is described by the PDE which I wrote down in this question. I am using an implicit numerical scheme to model the heat distribution evolving in time.

After some time has passed the system comes to a stationary distribution (corresponding with thermal equilibrium). One of the most interesting features of the model is the nonuniqueness of stationary state: if we'll start from another initial distribution, then we may come to another stationary state.

So, one of the task I need to solve is to find the set of stationary states for a given set of parameters.

I am solving it with following algorithm: I run the simulation until the heat distribution stops changing (difference between current and past distribution is lesser than given epsilon). So, after one stationary distribution is found, I increase and decrease laser power to find the next stationary state. By varying initial distributions for given set of parameters I compute a set of stationary distributions.

But, there are two problems. First, to compute the needed set with given accuracy I need to use small steps in space and time, small value of epsilon used to find stationary states and many steps of laser power variating. This is very computationally intensive.

The second problem is that it's easy to miss a distribution or count one as two, and I am not sure if there are other distributions which my algorithm cannot found.

So, what algorithms should I consider to use instead of just using implicit scheme and waiting until equilibrium is estabilished? There are two significant cases:

• there is one possible stationary state, and the task is just to find it with given accuracy and minimal amount of computation

• there is a set of possible stationary distributions and the task is to find them all.