Solving a system of polynomials in $N$ variables Suppose I am given some non-negative constants $(C_p)_{p=1, ..., l}$ and I would like to find an integer $N$ and vector $v \in R^N$ such that
$$
\sum_{i=1}^N (v_i)^p = C_p
$$
for $p=1, ..., l$.


*

*Can I find $v$ constructively?

*How small can I take $N$?

 A: It seems that for a generic right-hand side the $N = l$, but I cannot prove that.
I used to numerically solve a similar problem using the following algorithm I call continuous point insertion. The idea is to construct solution gradually, starting from the simpliest system
$$
v_1 = c_1
$$
and adding one more equation at a time.
Assume that we've already solved the system with $k$ equations, so we know $\mathbf{v}^{(k)} \in \mathbb R^k$ satisfying system
$$
\sum_{i=1}^k v_i = C_1\\
\sum_{i=1}^k v_i^2 = C_2\\
\vdots\\
\sum_{i=1}^k v_i^k = C_k\\
$$
Now, let's insert a new variable $v_{k+1}$ along with the $k+1$ equation parametrically
$$\begin{align}
\sum_{i=1}^k v_i + v_{k+1} &= C_1 + (1-\gamma) A_{k+1}\\
\sum_{i=1}^k v_i^2 + v_{k+1}^2  &= C_2 + (1-\gamma) A_{k+1}^2\\
\vdots\\
\sum_{i=1}^k v_i^k + v_{k+1}^k &= C_k + (1-\gamma) A_{k+1}^k\\
\sum_{i=1}^k v_i^{k+1} + v_{k+1}^{k+1} &= \gamma C_{k+1} + (1-\gamma) A_{k+1}^{k+1}+ (1 - \gamma)\sum_{i=1}^k (v_i^{(k)})^{k+1} \end{align}
$$
The $A_{k+1}$ is some arbitrary parameter, an initial guess for $v_{k+1}$.
When $\gamma = 0$ the system reduces to the previous system for which the solution is known
$$
\mathbf{v}^{(k+1)}\big|_{\gamma = 0} = 
\left(v_1^{(k)}, v_2^{(k)}, \dots, v_k^{(k)}, A_{k+1}\right).
$$
When $\gamma = 1$ is it the system we want to solve. Let $\mathbf{v}(\gamma)$ be the smooth (in $\gamma$) solution to the system above.
Denoting the system as $\mathbf{F}(\mathbf{v}(\gamma)) = \mathbf{g}(\gamma)$ we can obtain a system of ODE for $\mathbf{v}(\gamma)$:
$$
\mathbf v'(\gamma) = \left(\frac{\partial \mathbf F(\mathbf v(\gamma))}{\partial \mathbf v}\right)^{-1} \mathbf g'(\gamma)
$$
with initial condition
$$
\mathbf v(0) = \left(v_1^{(k)}, v_2^{(k)}, \dots, v_k^{(k)}, A_{k+1}\right)
$$
This technique is known as numerical continuation.
Note that this ODE was extremely stiff and ill-conditioned for my problem, so I had to use double-double and quad-double precision to integrate it for sizes about $n = 20$.
The method of integration I used was a combination of explicit Euler method with Newton refinement at each step:

$s = 0, \quad \gamma_s = 0, \quad h = h_\text{initial}$
while $\gamma_{s} < 1$ do

$\gamma_{s+1} = \gamma_s + h$
$\tilde{\mathbf{v}} = \mathbf{v}(\gamma_{s}) + \left(\dfrac{\partial \mathbf F(\mathbf v(\gamma_s))}{\partial \mathbf v}\right)^{-1} \mathbf g'(\gamma_s)$
Solve $\mathbf F(\mathbf v) = \mathbf g(\gamma_{s+1})$ for $\mathbf v$ using $\tilde {\mathbf{v}}$ as an initial guess.
If Newton method failed to solve the equation, decrease $h$ and restart the step from the beginning.
Use the solution as $\mathbf v_{s+1}$. Adjust $h$ according to the number of needed Newton iterations. $s = s + 1$.


When we reach $\gamma = 1$ we solved the system with $k+1$ equation. Repeat the step until the system of $l$ equations is solved.
Note that both the Euler step and Newton refinement use the same Jacobi matrix
$$
\frac{\partial F_p}{\partial v_i} = p v_i^{p-1}
$$
which is not singular unless some $v_i$ are the same (it's a variation of the Wandermonde matrix). So we need to avoid crossing of the $v_i(\gamma)$ trajectories. Careful choice of $A_{k+1}$ may help here.
