In order that the surface was more like your picture I did:

1) Removed the plane  x=0.

2) Changed the sign in front of  z -term.

3) Changed the scaling along the  y-axis.

3) Added the plotting of the Bifurcation Set (op(eliminate({f, implicitdiff(f, y, x)}, x)[2]) = 0  is the projecting surface onto x=-7):

f := x^3-x*z+(1/3)*y:

A := plots[implicitplot3d](f, x = -7 .. 5, y = -50 .. 50, z = -5 .. 7, axes = boxed, grid = [60, 60, 60], style = surface, color = khaki, orientation = [175, -20, 15], lightmodel = light1):

B := plots[intersectplot](x+7 = 0, op(eliminate({f, implicitdiff(f, y, x)}, x)[2]) = 0, x = -7 .. 5, y = -50 .. 50, z = -5 .. 7, color = red, thickness = 4, linestyle = solid):  #  Bifurcation Set

C := plots[implicitplot3d](x = -7, x = -7 .. 5, y = -50 .. 50, z = -5 .. 7, style = surface, color = "LightGrey"):  # The plane  x=-7

plots[display](A, B, C); 

 Edited.

Here's an example of another way to solve this problem. This universal method allows you  from a list of certain objects to select objects that satisfy certain properties. I believe that all beginners must master this method:

restart;

A:=[2,4,6,8,10,12,8,6]:

select(i->A[i]=6, [$ 1..nops(A)]);

                                   [3, 8]

Procedure  SM  solves your problem:

SM:=proc(n)

local L;

uses combinat;

L:=permute([(-1)$n,1$n],n);

interface(rtablesize=infinity);

Matrix(L);

end proc:

Example of use:

SM(4);

Period  procedure finds the least positive period for all expressions of the form a linear combination of  {sin(a*t+a0), cos(b*t+b0), exp((c*t+c0)*I)}  with rational coefficients  a, b, c  and  a0, b0, c0  are any real constants.

Period:=proc(Expr)

local Expr1, S, T, T1, d;

Expr1:=evalc(subs(op(indets(Expr,symbol))=t, Expr));

S:=indets(Expr1,'trig') minus indets(Expr1,'realcons');

T:=map(s->coeff(op(1,s),t), S);

T1:=map(t->2/t,T);

d:=ilcm(op(denom~(T1)));

ilcm(op(d*~(T1)))/d*Pi;

end proc:

Examples of use:

Period(sin(t/2+sqrt(2))-4*cos(2*t/3-5));

                            12*Pi

Period(2+sin(3*x)+cos(2*x)-sin(x/2)+exp(7/2*I*x));

                             4*Pi

Edited.

Estimation  N=62/eps  is universal, but very overpriced. A direct calculation gives us the exact value of  integer positive  N  for each epsilon.

restart;

x:=n->n^2/(n^2+31*n+228);

x0:=limit(x(n), n=infinity);

N:=epsilon->max(1, floor(op(1, solve(abs(x(n)-x0)<=epsilon, n)[-1]))+1):

Example:

N(10^(-3));

                                  30977

Check:

evalf(abs(x(30977)-x0));

evalf(abs(x(30976)-x0));

                       0.0009999791325

                        0.001000011390

Edited. 

To find stationary points use  fsolve  as in the simple example:

deriv:=x->fdiff(2*t-t^2, t=x);

fsolve('deriv'(x)=0);

                                  1.000000000

or the most simple

fsolve(deriv);

                                  1.000000000

Just write:

restart;

vec:=<a, b, c>:

a:=1:

vec;

 Addition: If you don't know the symbol for  vec[1]  then you can write:

restart;

vec:=<a,b,c>:

assign(vec[1], 1);

vec;

a;

Error reason is that you miss the multiplication sign between  b  and subsequent expression in parentheses. In this case, Maple treats  b  as the function name, and the expression in parentheses as the argument of this function. Here is another example where this situation is evident:

eval(sin(x), sin=1);

                                            1

If you put a multiplication sign, the result is true:

c:=1/(b*(1+x/l));

eval(c, [b=1,l=1]);

Another method with similar idea but slightly shorter:

eq := (-Omega^2*a*A[2]-Omega^2*m*B[1]+Omega*A[1]*c[1]+B[1]*k[1])*cos(Omega*t)+(Omega^2*a*B[2]-Omega^2*m*A[1]-Omega*B[1]*c[1]+A[1]*k[1])*sin(Omega*t) = 0:

eq1:=subs({sin(Omega*t)=s, cos(Omega*t)=c}, lhs(eq));

coeff(eq1, s);

coeff(eq1, c);

plots[inequal]  command allows you to get a high-quality picture:

plots[inequal]({y>=0, x^0.5+y^0.5<=2}, x=0..1, y=0..4, optionsfeasible=[color=yellow], optionsclosed=[color=blue, thickness=4], optionsexcluded=[color="WhiteSmoke"]);

 Edit. The above variant works in Maple 2015. For older versions Maple use  filledregions=true  option. The result is the same:

A:=plot([solve(x^0.5+y^0.5-2, y), 0], x=0..1, color=blue, thickness=4):

B:=plots[implicitplot](x^0.5+y^0.5-2, x=0..1, y=0..4, coloring=[yellow, "WhiteSmoke"], filledregions=true, gridrefine=3):

plots[display](A, B);

tickmarks  option solves your problem.

Example:

A:= LinearAlgebra[RandomMatrix](8, generator=5..15):

plots[matrixplot](A, axes=normal, tickmarks=[[seq(i=eq||i,i=1..8)], [seq(i=q||i,i=1..8)], default], view=[0..9, 0..9, 0..17], axesfont=[TIMES,ROMAN,10], orientation=[25,75]);

The symbol  r  is your fraction.

Coefficients_of_Fractions1.mw

applyrule([

cos(u::anything)*cos(v::anything)-sin(u::anything)*sin(v::anything)=cos(u+v),

cos(u::anything)*sin (v::anything)+sin(u::anything)*cos(v::anything)=sin(u+v),

sin(u::anything)*sin(v::anything)-cos(u::anything)*cos(v::anything)=-cos(u+v),

-sin(v::anything)*cos(u::anything)-sin(u::anything)*cos(v::anything)=-sin(u+v)], simplify(-(sin(psi[1](t))*cos(theta[1](t))*cos(gamma[1](t))+sin(psi[1](t))*sin(theta[1](t))*sin(gamma[1](t))-cos(theta[1](t))*cos(psi[1](t))*sin(gamma[1](t))+cos(psi[1](t))*sin(theta[1](t))*cos(gamma[1](t)))*cos(beta[1](t)), size)):

(trigsubs(2*combine(%))[])/2;

                         cos(beta[1](t))*sin(-theta[1](t)-psi[1](t)+gamma[1](t))

There are infinitely many solutions. They all lie on the circle with center  -I  and radius 1. The equation of the circle in the complex plane will:

x:=a+I*b:

evalc(abs(I+x) = 1);

                          (a^2+(b+1)^2)^(1/2) = 1

The real solution is only one  x=0

Addition: This result is clear without any calculations because the equation can be rewritten as  abs(x-(-I))=1  and we know that  abs(z1-z2)  is the distance  between  z1  and  z2  in the complex plane.

Perhaps you work in a classic interface, in which the filling for non-convex polygons is not working properly. In the standard worksheet everything correctly:

L := [[[0, 0, 0], [0, 1, 0]], [[x, x^2+1, 0], x = 0 .. 2], [[2, 5, 0], [2, 2, 0]], [[x, x, 0], x = 2 .. 0]]:

plots[display](Picture(L, color = yellow), axes = normal, scaling = constrained);