It seems Maple doesn't have a built-in procedure for constructing a vector field on a surface. However, if we know the parametric equations of the surface and the components of the vector field as functions of points on the surface, then using the  plots:-arrow command we can easily construct any number of field vectors.

In the example below, several vector field vectors are plotted on the unit sphere. In this example, all vectors are perpendicular to the surface and have a length of 1/4.

restart;
S:=plots:-display(plottools:-sphere(style=line, color=black, numpoints=400), axes=none):
L:=[sin(s)*cos(t),sin(s)*sin(t),cos(s)]:
F:=seq(seq(plots:-arrow(L, L/4, color=red), s=0..3.14,3.14/6), t=0..6.28,6.28/12):
plots:-display(S, F);

Answer questions a) - e) easily manually (without Maple). 
Below is an animation of the trajectories and velocity vectors of these yachts:

restart;
P:=proc(t)
local Tr1, Tr2, Ya1, Ya2, V1, V2, v;
uses plots, plottools;
Tr1:=plot([4+3*s,5-4*s,s=0..t], color=red, thickness=3);
Tr2:=plot([1+8*s,-8+6*s,s=0..t], color=blue, thickness=3);
Ya1:=disk([4+3*t,5-4*t],0.4, color=red);
Ya2:=disk([1+8*t,-8+6*t],0.4, color=blue);
V1:=arrow([4+3*t,5-4*t],[3,-4], width = [0.08, relative], head_length = [0.2, relative], color=red);
V2:=arrow([1+8*t,-8+6*t],[8,6],  width = [0.04, relative], head_length = [0.1, relative], color=blue);
v:=textplot([[4+3*t+3,5-4*t-4,v1,align=[right,above]],[1+8*t+8,-8+6*t+6,v2,align=[right]]], font=[times,bold,16]);
display(Tr1,Tr2,Ya1,Ya2,V1,V2,v);
end proc:

plots:-animate(P,[t], t=0..4, frames=40, scaling=constrained, size=[600,600]);

We can easily find the general solution to this problem (not just the smallest solution) if we use the  isolve  command. The notation and the beginning are the same as vv's ones.

restart;
sol:=isolve(10*u+7=(7*10^n+u)/3); # General solution for _Z1::nonenegint
eval(sol, _Z1=0);
eval(k=7*10^n+u, %); # Smallest solution

Since the value of each variable obviously does not exceed 637, the simplest enumeration in a double loop in a couple of seconds gives all 8 solutions:

restart;
for x from 0 to 637 do
for y from 0 to 637 do
if sqrt(x)+sqrt(y)=sqrt(637) then print({'x'=x,'y'=y}) fi;
od: od: 

The last line should be:

plots[odeplot](soln, [[t, r(t)], [t, varphi(t)]]);

The task is not difficult and comes down to solving a simple equation:

restart;
# A is the center of one of the circles k3
A:=<1+r,0>: B:=<sqrt(2)*cos(Pi/5),sqrt(2)*sin(Pi/5)>:
# The distance between the points A and B must be equal r
sqrt((B-A).(B-A))=r assuming real; 
r:=solve(%);
rationalize(r);
expand(%);
evalf(%);

# Visualization

with(plots): with(plottools):
k1:=circle(color=blue): k2:=circle(sqrt(2), color=green):
k0:=circle([1+r,0],r, color=red):
k3:=seq(rotate(k0,2*Pi*i/5), i=0..4):
display(k1,k2,k3);

The reason Maple doesn't fully expand the parentheses is because you forgot to put a multiplication sign between the parentheses in the middle term. If you want Maple to write the resulting two-variable  с  and  m  polynomial as a one-variable  m  polynomial, you should use the  sort  command:

If you omit this multiplication sign, as you did, Maple interprets this entry in the middle term as calling the sum of two functions  5/2 (constant function)  and 12 c  applied to the argument m-1/2, followed by squaring.

Another example of a function call:

(a+3*b+4)(c+d);

You forgot to put the multiplication sign between the brackets in the middle term.

I was about to answer, but Preben Alsholm beat me to it.

In my opinion, Maple's behavior is quite logical. If you do not specify specific initial conditions, then Maple formally returns 2 branches of solutions. But in each specific case, only one branch should be chosen depending on the sign of  v__0. Since you do not specify the meaning of  v__0, then Maple simply does not know which branch should be returned.

restart;	
ode:=v(x)*diff(v(x),x) = g;
dsolve(ode);
ic:=v(x__0) = v__0;
sol1:=dsolve([ode,ic]) assuming v__0>0;
sol2:=dsolve([ode,ic]) assuming v__0<0;

Download A_way.mw

I think it would be more interesting to animate the curve  y = sin(x)^n  shape changing as the exponent  n  changes. Below in the code the exponent increases from  1/6  to  6 :

restart;
with(plots):
L:=[surd(sin(x), 6), surd(sin(x), 5), surd(sin(x), 4), surd(sin(x), 3), surd(sin(x), 2), sin(x), sin(x)^2, sin(x)^3, sin(x)^4, sin(x)^5, sin(x)^6]:
T:=[seq(textplot([3*Pi/8,eval(t,x=3*Pi/8)+0.06,y=t], font=[times,bold,16], align=left), t=L)]:
L1:=[seq(display([plot(t[1],x=0..Pi, color=red, thickness=3, scaling=constrained),t[2]])$10, t=zip((x,y)->[x,y],L,T))]:
display(L1, insequence, size=[800,400]);

I assumed that  i  (in Maple it should be coded as I )  is the imaginary unit, and  w  is a real parameter.

restart;
Expr:=(2500*I*w)/(1+5*I*w)+(200*I*w)/(1-10*I*w)+5;
evalc(Expr);

# Or

normal(Re(Expr))+normal(Im(Expr))*I assuming real;

restart;
kernelopts(version);

ode:=diff(y(x),x) = (2*sin(2*x)-tan(y(x)))/x/sec(y(x))^2:
ode1:=simplify(ode);

dsolve(ode1);
dsolve(ode1,implicit);

To avoid indexing issues (for matrices, indices start at 1, not 0), we can use arrays instead of matrices. See the visualization below for an example from the   wiki  https://en.wikipedia.org/wiki/Multinomial_distribution  :

restart;
M:=Array(0..6,0..6,{seq(seq((i,j)=6!/(i!*j!*(6-i-j)!)*0.2^i*0.3^j*0.5^(6-i-j),
j=0..6-i),i=0..6)}):
plots:-matrixplot(M, heights = histogram, gap=0.5, tickmarks=[[seq(i=i-1/4,
i=1/4..25/4,1)],[seq(i=i-1/4,i=1/4..25/4,1)],default], labels=[i,j,P(i,j)], 
labelfont=[times,16]);

Use the seq command in plots:-textplot .

For some reason the worksheet doesn't open here - so I'm only providing a link to download it.

Download S7MAA_Dveloppement_Limit_new.mw

The final picture: