Solve system for diagonal matrix In trying to write an integral relation in a discrete manner, I got to an equation of the form
$$MAM x=b$$
where $A$ is a given symmetric  matrix, $b$ and $x$ are given vectors and unknown $M$ is a diagonal matrix. The values of all the matrices and vectors can be complex. How can I solve this for $M$?

My attempts
I tried to find patterns that would allow for inverting the order in which the matrices are applied but without a solution. Another thing I tried has been to actually write down explicitly the relation in order to get the system of equations. For the first equation, I got something like
$$m_{11} \left( a_{11}m_{11}x_1 + a_{12}m_{22}x_2 + \cdots + a_{1j}m_{jj}x_j + \cdots + a_{1n}m_{nn} x_n \right) = b_1$$
where $m_{ii}$ are the diagonal entries of $M$, $a_{ij}$ are entries of $A$ and so on for vectors $x$ and $b$, and so on for the other equations. This system of equations is not linear and I have no idea how to solve a non-linear system of equations.
I assume that another way in which the problem can be formulated, although it is not really the same is how can one solve a system of equations like
$$Ax = b/x$$
where $b/x$ should be interpreted as a vector defined by the element-wise division of each pair of points $(b/x)_j = b_j/x_j$.
 A: $
\def\o{{\tt1}}\def\p{\partial}
\def\L{\left}\def\R{\right}
\def\LR#1{\L(#1\R)}
\def\BR#1{\Big(#1\Big)}
\def\BBR#1{\Bigg(#1\Bigg)}
\def\vec#1{\operatorname{vec}\LR{#1}}
\def\diag#1{\operatorname{diag}\LR{#1}}
\def\Diag#1{\operatorname{Diag}\LR{#1}}
\def\trace#1{\operatorname{Tr}\LR{#1}}
\def\qiq{\quad\implies\quad}
\def\grad#1#2{\frac{\p #1}{\p #2}}
\def\dk#1{\LR{#1_k-#1_{k-\o}}}
\def\c#1{\color{red}{#1}}
$Let's rename the variable $M\to Y,\,$ then the independent variable for this problem is the vector $y$
$$\eqalign{
Y &= \Diag{y} \quad\iff\quad y=\diag{Y} \\
}$$
Define the vector-valued function
$$\eqalign{
f(y) &= \BR{YAYx-b} 
  &= \BR{A\odot{yy^T}}x - b \\
}$$
where $(\odot)$ denotes the elementwise/Hadamard product.
The current problem can be restated as
a system of (mildly) nonlinear equations
$${f(y) = 0}$$
The simplest route to a solution is probably the Barzilai-Borwein method, since (in this situation) it does not require any gradient calculations.
The basic Barzilia-Borwein method is extremely straightforward
Initialize
$$\eqalign{
y_0 &= random
 \qquad\qquad\qquad\qquad\qquad\quad \\
}$$
First step
$$\eqalign{
f_0 &= f(y_0) \\
y_1 &= y_0 - \LR{\frac{0.0\o\,\|y_0\|}{\|f_0\|}}f_0
 \qquad\qquad\quad\quad \\
k &= \o \\
}$$
Subsequent steps
$$\eqalign{
f_k &= f(y_k) \\
y_{k+\o} &= y_k - \LR{\frac{\dk{y}^T\dk{f}}{\dk{f}^T\dk{f}}}f_k \\
k &= k+\o \\\\
}$$
Stop when the function is nearly zero $\;\|f_k\|\le 10^{-12} $
or when the steplength gets too small
$\;\BBR{\frac{\|y_k-y_{k-\o}\|}{\o+\|y_k\|}}\le 10^{-12}$


The introduction of a per iteration difference operator
$$\eqalign{
dy_k &= y_k - y_{k-\o} \\
}$$
allows for a more concise description of the algorithm
$$\eqalign{
y_{k+\o} &= y_{k} - \LR{\frac{dy_{k}^Tdf_{k}}{df_{k}^Tdf_{k}}}f_{k}, \qquad
f_{k+\o} &= f(y_{k+\o}) \\
}$$
Suppressing the $k$-index makes it look less cluttered
$$\boxed{\eqalign{
y_+ &= y - \LR{\frac{dy^Tdf}{df^Tdf}}f,
 \qquad
 \qquad
f_+ &= f(y_+) \\
}}$$
