Dear Debika,

While I am not a particular expert on the specific problem you are dealing with, to my knowledge self-consistent equations that arise in quantum mechanics are typically solved by using a fixed-point iteration (heuristically speaking at least). If x is the object you wish to solve, you formulate the self-consistent problem in the form x = F(x) and then solve the iteration xn+1 = F(xn). Here F could be, e.g., a Hamiltonian that has some terms whose evaluation requires you to solve a (typically linear) subproblem for a fixed xn. You start with an initial guess x0, then solve the context-dependent subproblem x1 = F(x0) where you assume x0 to be known, and repeat. By continuing in this way, the iteration hopefully converges to a fixed point, such that the norm |xn+1-xn| has becomes smaller than your tolerance for error.

Codingwise, adding this "topmost fixed-point iteration layer", which is there to enforce the self-consistency, is straightforward if you have the solver for the context-dependent subproblem. Depending on the number of iterations you need, it can of course increase the need for computation time significantly.

If you are dealing with self-consistent equations in terms of Wigner functions, then you could try using spectral split propagator, which is a very efficient method and there are python programs written for it, e.g. here: http://www.quantuminfodynamics.com/ (see section "Software").

Tigran Aivazian Thank you for your reply. My equation doesn't contain Wigner function. Therefore, I am still stick with the problem!

Pekko Kuopanportti Thanks a lot for your description in details.