I don't know what the azimuth and polar angles are.  In the attached worksheet I have used the longitude and colatitude angles which are normally used in spherical coordinates.  Perhaps they mean the same thing.

spherical-cap.mw

The problem with your plot command is that it does not refer to sol which contains the numerical solution of the differential equation. What you want is
   
or equivalently
   

You may add time range and color options, as needed.

One solves an equation.  One evaluates an expression.  I think you mean evaluate, not solve.

Anyway, here it is:

evalf(BesselI(0, 0.012*(-1)^(1/4)*sqrt(1000)));
        

I am guessing that you encoded ϕ_c as ϕ[c]. That's the source of your trouble. Instead, encode it as ϕ__c, that is, ϕ, followed by two underscores, followed by c.

Since you have not specified the parameter values, I will take your inequality to be:

and let q = 3.

restart;
eq := 1 + 2*q - 2*K*(1+q)/N;
q := 3;
tmp := plot3d(eq, K=0..1, N=0..1, view=0..8,
   axes=boxed, labels=[K, N, Gamma],
   filled=[style=surface, transparency=0.2]);
plots:-display(tmp, view=0.01..8);

Your worksheet is presented as a horrible mess.  I did not not have the courage to wade through it completely.  It appears that it is suffering from the 2Dmath-hormone-imbalance-syndrome.  I suggest that you consider switching to 1D input, as described in Configureing Maple, although this may be too late for you since you need to start from a fresh worksheet—I don't know how to rescue an existing worksheet.

Anyway, here is the substitution that you are looking for:

subs(
    x__s_dot(t)=diff(x__s(t),t),
    y__s_dot(t)=diff(y__s(t),t),
    ...
    u__2_dot(t)=diff(u__2(t),t),
  xEQNS):

I have shown only three of the fifteen substitutions here.  You may add the missing twelve yourself.

As we see in mehdi jafari's response, the exact solution is too complicated to be useful.  A numeric solition is more useful in this case, but for that you need to specify initial conditions, such as:

ic := A1(0)=1, D(A1)(0)=2, A2(0)=3, D(A2)(0)=4, A3(0)=5, D(A3)(0)=6;
sol := dsolve({eq1, eq2, eq3, ic}, numeric);
plots:-odeplot(sol, [[t, A1(t)],[t, A2(t)],[t, A3(t)]], t=0..6, color=[red, green, blue]);
 

By using textplot(), you may place a text anywhere you want.  For example:

restart;
with(plots):
myfont := [Times,24]:
p1 := plot(x^2, x=-1..1);
p2 := textplot([-0.2, 0.8, typeset(y=x^2)], align=left, font=myfont);
display([p1,p2], labelfont=myfont);

I can't tell what you did to get to what you have shown.  Copy and paste the following into Maple—it should do what you want:
PDE := diff(u(x,t),t) = diff(u(x,t),x,x) + 2*u(x,t)^2*(1 - u(x,t));

It is not difficult to show that the general solution of your differential equation is

where A and B are arbitrary constants.  For continuity at x=1 you need A=0, but then that leaves you with only B to play with, which means that you cannot specify two boundary conditions idependently.

Try these:

odeplot(sol, [t, Re(z(t))], t = 0 .. 3);
odeplot(sol, [t, Im(z(t))], t = 0 .. 3);
odeplot(sol, [Re(z(t)), Im(z(t))], t = 0 .. 3);

The proper way of integrating a function over a triangle is through the triangle's barycentric coordinates.  The worksheet integration-over-triangle.mw sets up a general framework for integrating any function f(x,y) over an arbitrary triangle.  At the bottom of the worksheet the general framework is applied to your specific problem.

Note added on August 22, 2017:

Robert Lopez has pointed out to me that the worksheet's computation sometimes gets the sign of the integral wrong.  This is because the worksheet assumes implicitly that the triangle's vertices are enumerated in the counter-clockwise order.  To get the correct result regardless of the ordering of the vertices, replace A with abs(A) in the expression for T_integral.

The expression which you wish to expand is a fraction whose numerator and denominator vanish with high order at [u,v]=[0,0].  A good strategy in this case is to expand the numerator and denominator separately, then combine the results, like this:

restart;
expr := ... your messy expression here ...;
Num := numer(expr):
Den := denom(expr):
N := 19;    # order of expansion
Num_N := mtaylor(Num, [u,v], N);
Den_N := mtaylor(Den, [u,v], N);
mtaylor(Num_N/Den_N, [u,v], N):  collect(%, [u,v], distributed);

The result is:

Increase N to get more terms.

Here is a self-contained worksheet for your experimentation: worksheet.mw

In your for-loop you have:

The values of all variables on the right-hand side are specified before you enter the for-loop.  Why do you expect them to change as you itererate through the loop?

Since the right-hand side is fixed, all x₁₁ have the same value, and that's why your point does not move.  You need to re-think the strategy.

An equation of the type L(n) = a*L(n-1) + b*L(n-2) is called a recurrence relationship of order two or a difference equation of order two.  To be overly pedantic, it would be called a linear homogeneous difference equation of order two with constant coefficients.

Maple's rsolve() can solve such equations explicitly.  For example:

eq := L(n) = a*L(n-1) + b*L(n-2);
sol := rsolve({eq, L(1)=p, L(2)=q}, L(n));

In your worksheet you have a=1, b=10, p=6, q=2.  We obtain the corresponding solution through:
mysol := eval(sol, {a=1, b=10, p=6, q=2});

Here we produce the first 10 terms of the sequence:
simplify([seq(mysol, n=1..10)]);
    [6, 2, 62, 82, 702, 1522, 8542, 23762, 109182, 346802]