Thank you, Jean-Mark, now everything is ok. I have only a little question: when the program needs a very long time, how can I know (after some hours) how many hours the computer needs more, to finish the calculus?

Thanks, gulliet, for the correct paranthese. But the big problem is that I can not plot the curves. One curve for one value of alpha. I need to plot 157 curves, for every value of paramether alpha. Can you help me, please? The program now is:

restart:with(plots):
Warning, the name changecoords has been redefined

> mu:=0.1:
> rs:=((x+mu)^2+y^2)^0.5:
> re:=((x+mu-1)^2+y^2)^0.5:
> Q:=beta/10*(1-mu)*(cos(alpha/100))^2:
> S:=(x+mu)*cos(alpha/100)+y*sin(alpha/100):
> T:=-(x+mu)*sin(alpha/100)+y*cos(alpha/100):
> eq1:=(1-mu)*(x+mu)/rs^3+mu*(x-mu+1)/re^3-x=Q*S/rs^3:
> eq2:=(1-mu)*y/rs^3+mu*y/re^3-y=Q*T/rs^3:
> S :=[seq([seq(fsolve({eq1,eq2},{x,y}), beta=0..10)],alpha=0..157)]:
> S[1]:
> subs(S[1], [x,y]):
Error, wrong number (or type) of parameters in function subs

> p:=seq(subs(i, [x,y]), i=S):
> p:=pointplot({p},axes=boxed,symbol=circle, color=black):

I use now the ideea of roman_pearce to plot the solution of my system. It is ok, thanks roman_pearce. The parameter was beta and the system was solved for a very particular value of parameter alpha. I obtained some points, situated on a curve. But now I want to give 100 values to parameter alpha and obtain points situated on 100 curves. I tried to convince Maple to do the job, but ...ERROR. Can somebody help me, please? The program is :

restart:with(plots):
Warning, the name changecoords has been redefined

> mu:=0.1:
> rs:=((x+mu)^2+y^2)^0.5:
> re:=((x+mu-1)^2+y^2)^0.5:
> Q:=beta/10*(1-mu)*(cos(alpha/100))^2:
> S:=(x+mu)*cos(alpha/100)+y*sin(alpha/100):
> T:=-(x+mu)*sin(alpha/100)+y*cos(alpha/100):
> eq1:=(1-mu)*(x+mu)/rs^3+mu*(x-mu+1)/re^3-x=Q*S/rs^3:
> eq2:=(1-mu)*y/rs^3+mu*y/re^3-y=Q*T/rs^3:
> S :=[seq([seq(fsolve({eq1,eq2},{x,y}), beta=0..10))],alpha=0..157]:
Error, `)` unexpected

> S[1]:
> subs(S[1], [x,y]):
> p:=seq(subs(i, [x,y]), i=S):
> p:=pointplot({p},axes=boxed,symbol=circle, color=black):

Thanks for helping me. I used the solution of Ken.Lin @Taiwan, but with the value of parameter mu changed, the eval instruction create me troubles. Please , what is wrong now?

The program is now:

restart:with(plots):
> alpha:=Pi/3:
> mu:=3.006*10^(-6):
> rs:=((x+mu)^2+y^2)^0.5:
> re:=((x+mu-1)^2+y^2)^0.5:
> Q:=beta/10*(1-mu)*(cos(alpha))^2:
> S:=(x+mu)*cos(alpha)+y*sin(alpha):
> T:=-(x+mu)*sin(alpha)+y*cos(alpha):
> eq1:=(1-mu)*(x+mu)/rs^3+mu*(x-mu+1)/re^3-x=Q*S/rs^3:
> eq2:=(1-mu)*y/rs^3+mu*y/re^3-y=Q*T/rs^3:
> plist:=[]:
>    for beta from 0 to 10 do
>       fsolve({eq1,eq2}, {x,y});
>       xx:=eval(x,%):
>       yy:=eval(y,%%):
>       plist:=[op(plist),[xx,yy,beta/10]]:
>    end do:
Error, wrong number (or type) of parameters in function evalf

> plist;
> plots[pointplot3d](plist,axes=boxed,symbol=point, color=red):