Easier way to solve equation systems of $a+b+c+\cdots{}= 1$, $a^2 + b^2 + c^2+\cdots{}=2$ and so on without having to crunch massive expressions I study at below college level. I have been trying to solve certain systems of equations involving $n$ equations of $n$ unknowns. For example, for $2$ unknowns, the problem is
\begin{align} 
a^{\phantom{1}} + b^{\phantom{1}} &= 1 \\ 
a^2 + b^2 &= 2 \\ 
a^3 + b^3 &={} ?
\end{align}
This can be solved with elementary algebra and/or WolframAlpha. You can generalize this to more unknowns:
\begin{align} 
a^{\phantom{1}} + b^{\phantom{1}} + c^{\phantom{1}} &= 1 \\ 
a^2 + b^2 + c^2 &= 2 \\ 
a^3 + b^3 + c^3 &= 3 \\ 
a^4 + b^4 + c^4 &={} ? 
\end{align}
with the same restraint: $n$ unknowns, $n$ equations, in each equation the powers of each variable is the same, and the pattern is clear.
Now, I, off of only the first $3$ cases (including the trivial case $a = 1$, find $a^2$) made a conjecture about the result (the missing value of the final expression). Since this is such a random guess at the value, and so many functions could meet just the first few data points, I want to solve the version with $4$ unknowns, just to see whether the conjecture's true.
However, this is very difficult. The expansions quickly get out of hand and not even WolframAlpha can do it. I want a way to at least get the solving process under control. Usually, one'd generate equations and use those to solve for things like $a \cdot b^3$, but here the issue is that just setting up the equations is a task too difficult. Is there a way to elegantly solve the system? I don't mind trading in time for maybe some more difficult math.
 A: Actually, Mathematica can solve your equations. The
code
x/.Solve[{a+b+c==1, a^2+b^2+c^2==2, a^3+b^3+c^3==3,
     a^4+b^4+c^4==x},{a,b,c,x}] //InputForm

returns the result $\,x=25/6\,$ with
{25/6, 25/6, 25/6, 25/6, 25/6, 25/6}

The key idea for an elegant solution is to use Newton's identities. As Wikipedia states

In mathematics, Newton's identities, also known as the Girard–Newton formulae, give relations between two types of symmetric polynomials, namely between power sums and elementary symmetric polynomials.

Thus, we are given the power symmetric polynomials
$$ p_1 = a+b+c = 1,\\ p_2 = a^2+b^2+c^2 = 2,\\
p_3 = a^3+b^3+c^3 = 3. $$
Use Newton's identities to find the elementary symmetric polynomials to get
$$ e_1 = p_1 = 1, \\
e_2 = (p_1^2-p_2)/2 = -\frac12, \\
e_3 = \left(\frac12 p_1^2 - \frac32 p_1p_2 + p_3\right)/3 = \frac16. $$
There is a recursion for the $\,p_i\,$ given by
$$ p_{i+3} = p_{i+2}e_1 -p_{i+1}e_2 +p_{i}e_3. $$
Apply this for $\,i=1\,$ gives the result
$$ p_4 = p_3\cdot(1) -p_2\cdot(-1/2) +p_1\cdot(1/6) = \frac{25}6. $$

You need to use the relation between the
first $\,n\,$ elementary symmetric functions and the
first $\,n\,$ power symmetric functions to solve for
$\,p_{n+1}\,$ given $\,p_1,p_2,\dots,p_n.\,$ All of
the formulas generalize.
However, in the special case were $\,p_i=i\,$ for
$\,i=1,\dots,n\,$ the answer is already given in the
OEIS sequence A094094

Define x[1]...x[n] by the equations Sum_{j=1..n} x[j]^i = i, i=1..n; a(n) = n! * Sum_{j=1..n} x[j]^(n+1).

The entry gives an exponential generating function

E.g.f.: (1-exp(x/(x-1)))/(1-x)^2.

and also closed form formulas and summations

a(n) = n!*(n+1-LaguerreL(n,1,1)) = Sum_{k=1..n} (-1)^(k+1)*n!/k!*binomial(n+1,k+1).

