Well, you did get to a start, and this does not look directly like homework.
I'll rewrite your function, making it a function of only S1, then call solve inside.
function res = myfun(S1)
syms x1 x2 x3 x4
w1 = 1;
w2 = 3;
a1 = -24.851;
a3 = -66.319;
a2 = -9.247;
a8 = -345.019;
a4 = 0;
a6 = -66.319;
a5 = -3.082;
a7 = 0;
mu1 = 0;
mu2 = 0;
f1 = 0;
f2 = 0;
S2 = 0.01;
E1 = 8 * w1 * mu1 * x1 - a2 * x1^2 * x2 * sin(x3) - 4 * f1 * sin(x4);
E2 = 8 * w2 * mu2 * x2 + a5 * x1^3 * sin(x3);
E3 = 8 * w1 * x1 * S2 + (3* a1 * x1^2 + 2 * a3 * x2^2) * x1 + a2 * x1^2 * x2 * x2 * cos(x3) + 4* f1 * cos(x4);
E4 = 8 * w2 * x2 * (3 * S2 - S1) + (3 * a8 * x2^2 + 2 * a6 * x1^2) * x2 + a5 * x1^3 * cos(x3);
res = solve(E1,E2,E3,E4,x1,x2,x3,x4);
It looks like there are generally 9 independent solutions, 3 of them are not a function of S1 at all. As well, x(2) always seems to come out as zero in all solutions.
Lets take a look at a few cases for different values of S1.
res = myfun(-50);
vpa([res.x1,res.x2,res.x3,res.x4],5)
ans =
[ 0, 0, 0, -1.0771]
[ 0, 0, 0, 1.0771]
[ 0, 0, -0.032758, -9.0236e-8]
[ 0, 0, 0, 0]
[ 3.1416, 0, 0, -1.0771]
[ 3.1416, 0, 0, 1.0771]
[ 3.1416, 0, 0, 0]
[ 0, 0, 1.5708, 0]
[ 3.1416, 0, 0.032758, -9.0236e-8]
res = myfun(-30);
vpa([res.x1,res.x2,res.x3,res.x4],5)
ans =
[ 0, 0, 0, -0.83445]
[ 0, 0, 0, 0.83445]
[ 0, 0, 0, 0]
[ 3.1416, 0, 0.032758, -1.5035e-7]
[ 3.1416, 0, 0, -0.83445]
[ 3.1416, 0, 0, 0.83445]
[ 0, 0, -0.032758, -1.5035e-7]
[ 3.1416, 0, 0, 0]
[ 0, 0, 1.5708, 0]
res = myfun(0);
vpa([res.x1,res.x2,res.x3,res.x4],5)
ans =
[ 0, 0, 0, -0.026374]
[ 0, 0, 0, 0.026374]
[ 3.1416, 0, 0.032757, -0.00018753]
[ 0, 0, 0, 0]
[ 0, 0, -0.032757, -0.00018753]
[ 3.1416, 0, 0, -0.026374]
[ 3.1416, 0, 0, 0.026374]
[ 3.1416, 0, 0, 0]
[ 0, 0, 1.5708, 0]
res = myfun(25);
vpa([res.x1,res.x2,res.x3,res.x4],5)
ans =
[ 0, 0, 0.032758, -1.8073e-7]
[ 0, 0, 0, -0.76091i]
[ 0, 0, 0, 0.76091i]
[ 0, 0, 0, 0]
[ 3.1416, 0, 0, 0]
[ 0, 0, 1.5708, 0]
[ 3.1416, 0, -0.032758, -1.8073e-7]
[ 3.1416, 0, 0, -0.76091i]
[ 3.1416, 0, 0, 0.76091i]
res = myfun(10);
vpa([res.x1,res.x2,res.x3,res.x4],5)
ans =
[ 0, 0, 0, -0.48081i]
[ 0, 0, 0, 0.48081i]
[ 0, 0, 0, 0]
[ 3.1416, 0, 0, 0]
[ 0, 0, 1.5708, 0]
[ 3.1416, 0, -0.032758, -4.5249e-7]
[ 0, 0, 0.032758, -4.5249e-7]
[ 3.1416, 0, 0, -0.48081i]
[ 3.1416, 0, 0, 0.48081i]
Be careful, because many of the solutions are the same, just reported in a different sequence.
So x2 is always zero for all the cases I've tried. You could probably show that to be necessary by looking at the equations.
Those 3.1416 values are pi, rounded off to 5 digits by vpa. In fact, the symbolic toolbox shows it to be exactly pi. And pi/2 is 1.5708.
So some solutions in common to all are [0,0,0,0], [0,0,pi/2,0], and [pi,0,0,0]. I'd want to ignore all of those apparently trivial solutions.
Next, while all solutions seem to be real for negative S1, at some point it seems that solutions start turning imaginary for at least some values of positive S1.
There appear to be 4 solutions of the form
[0,0,0, x(4)]
[0,0,0,-x(4)]
[pi,0,0, x(4)]
[pi,0,0,-x(4)]
There x(4) is a function of S1. With a little algebra, that relationship could surely be uncovered.
Finally, there always seem to be a pair of solutions of the form
[ 0,0,-x3,x4]
[pi,0, x3,x4]
Here x3 and x4 are both functions of S1. Again, some simple algebra should recover that relationship. That seems to account for all 9 basic solutions.
Quick inspection of your coefficients has many of them as zero, f1, for example. A result of that is Equation 2:
E2
E2 =
-(1541*x1^3*sin(x3))/500
E2 is zero ONLY when x1 is zero, OR sin(x3) is zero. That happens at x1==0, x3==0, or x3==k*pi, for any integer k.
You can do the rest.
