# How to solve self-consistent equations numerically?

I am working in condensed matter, where I'm required to solve an integral self consistently which is of the form, $$\Delta = \int f(x,y,\Delta,\mu,h)dxdy$$ Basically I need to find value of $$\mu,\,\Delta$$ such that the above equation is satisfied for the specific value of $$h$$. Also the solution ($$\mu^*, \Delta^*$$) should be such that $$\int g(x,y,\Delta^*,\mu^*,h) dx dy=1$$ And the $$h$$ is a free parameter. How to solve such a problem numerically?

What I have tried until now is to take a range of $$\mu$$ and $$\Delta$$ and check if both the equations are satisfied. This is usually taking very very long time numerically. If anyhow I reach my solution how to know if that is the true minima of the solution space?

It seems to me that you're looking for constrained optimization. These are often reasonably similar to more common optimization methods, such as Nelder-Mead or Newton's method, except they usually add a penalty function that increases the 'cost' of a point (in $$\{\eta,\,\Delta\}$$ space) when the constraint is not satisfied (in your case, when the integral over $$g$$ is not 1).
We can treat this as a problem of solving 2 equations of 2 variables. If we define $$F_1 \left(\mu,\Delta\right)=\Delta-\int f\left(x,y,\Delta,\mu,h\right)dxdy$$ and $$F_2 \left(\mu,\Delta\right)=1-\int g\left(x,y,\Delta,\mu\right)dxdy$$ then the problem becomes finding a $$\left(\mu^{*},\Delta^{*}\right)$$ which is a simultaneous root of both $$F_1$$ and $$F_2$$. For some vector function $$\vec{F}:\mathbb{R}^{n}\to\mathbb{R}^{n}$$, we can find a root (i.e. solve $$\vec{F}\left(\vec{x}\right)=\vec{0}$$) with the multivariable Newton-Raphson method. Starting with some guess $$\vec{x}_{old}$$, calculate the Jacobian matrix $$J(\vec{x}_{old})$$ of $$\vec{F}$$. If an analytic expression for the Jacobian is known, we use that (evaluate it at $$\vec{x}_{old}$$). Otherwise, we compute the derivatives in the Jacobian as finite differences, $$\frac{\partial F_{i}}{\partial x_{j}}\left(\vec{x}_{old}\right)=\frac{F_{i}\left(\vec{x}_{old}+h\hat{e}_{j}\right)-F_{i}\left(\vec{x}_{old}\right)}{h}$$ (here $$\hat{e}_j$$ is the $$j$$ unit vector). Now we look for a $$\vec{\delta}$$ such that $$\vec{x}_{new}:=\vec{x}_{old}+\vec{\delta}$$ will be closer to a root. So we stipulate $$\vec{F}(\vec{x}_{new})=\vec{0}$$. Applying a linear approximation around $$\vec{x}_{old}$$, we have $$\vec{0}=\vec{F}(\vec{x}_{new})\approx \vec{F}\left(\vec{x}_{\text{old}}\right)+J\left(\vec{x}_{\text{old}}\right)\vec{\delta}$$, and hence to find $$\vec{\delta}$$ we must solve the matrix equation $$J\left(\vec{x}_{\text{old}}\right)\vec{\delta}=-\vec{F}\left(\vec{x}_{\text{old}}\right)$$. This is done with Gaussian elimination. Now $$\vec{x}_{new}=\vec{x}_{old}+\vec{\delta}$$ is our new guess, and we repeat this process until convergence.