# Finding all solutions of a multidimensional nonlinear equation

I have a problem in which I have two (though there can be more) multivariate Gaussian functions very close to each other and I want to find whether being next together causes the sum of these two Gaussian functions to be unimodal or bimodal (more than one peak). I have tried to solve this problem analytically so far; since the type of multivariate Gaussians I am using looks something like this:

$$f(x)=\frac{1}{\sqrt{(2\pi)^k|\Sigma|}}\exp\left(-\frac{1}{2}(x-m)^T\Sigma^{-1}(x-m)\right)$$

where x is a k dimensional vector and sigma is the k x k covariance matrix. The obvious thing for me to do was the take the derivative since I am trying to find the peaks, which gives me something like this: $$\frac{\partial f(x)}{\partial x}=f(x)\Sigma^{-1}(x-m)$$

Now I'm at a loss since I need to solve this equation for all the zeros and I don't know how to do that. I know that there are optimising solutions out there to find the solutions and since my initial guess would be near m so that could work well except I am trying to figure out if the derivative has more than one peak and where it is - that method only tells me where one is; therefore, I need a method to find all zeros/roots. Does anyone know how to go about solving either the initial problem statement or finding the roots of the derivative. (I work in Python, but since the question is so open ended any solution would be appreciated)