You an do this using the  plots:-textplot  command:
 

restart;
ff:=x->x^2:
A:=plot(ff(x), x = 0 .. 10, color = black, labelfont = [TIMES, 19], thickness = 1, size = [1200, 600], axes = boxed, linestyle = 5, labels = ['x', ``], axis = [gridlines = [20, color = black]]):
B:=plots:-textplot([0,104.9,"y  "], align=left, font=[times,19]):
plots:-display(A,B);

Output:  

It is advisable to solve your new example using the  geometry  package:

restart;
with(geometry):
point(C,0,0):
point(B,7,0):
solve({x^2+y^2=25,(7-x)^2+y^2=9,y>0}, explicit);
point(A,eval([x,y],%)):
triangle(T,[A,B,C]);
incircle(inc,T);
v:=coordinates(center(inc));
r:=radius(inc);
map(coordinates, [A,B,C]);
map(t->t-v, %);
plots:-display(plot([%[],%[1]], scaling=constrained), plot([r*cos(t),r*sin(t),t=0..2*Pi]));

This example is easily solved by direct calculation (without any equations), if we use the Heron formula and the formula for the radius of the circumcircle (I chose the same position of the triangle as Carl did). For the second example 3-4-5, everything is the same:

restart;
a,b,c:=3,5,7: p:=(a+b+c)/2:
S:=sqrt(p*(p-a)*(p-b)*(p-c)); R:=a*b*c/4/S;
cosA:=b/2/R;  cosB:=c/2/R;
# So we have the final answer
A:=[b*cosA-R,b*sqrt(1-cosA^2)];
B:=[c*cosB-R,c*sqrt(1-cosB^2)];
C:=[-R,0];

Any for loop is not needed here. As an animation parameter, I took the distance from the origin to the lower end of the ladder:

restart;
with(plottools):
with(plots):
co := blue:
animate(display,[line([0,sqrt(400-a^2)],[a,0], color = co, thickness = 3)], a=0..20, frames=60);

Addition - the usage of the  trace  option:

restart;
with(plottools):
with(plots):
co := blue:
animate(display,[line([0,sqrt(400-a^2)],[a,0], color = co, thickness = 3)], a=0..20, frames=60, trace=[10,20,30,40,50,57]);

The method below is slightly longer, but I think it will be clearer for OP:

a:=1/(i*sqrt(i+1)+(i+1)*sqrt(i));
b:=expand(rationalize(a));
c:=normal(`+`(op(1..2,b)))+op(3,b);
limit(sum(c, i = 1 .. n), n = infinity);

restart;
randomize():
A:=[1,1.732,1.23,4.42,9,6.45,3.45,8.428,9.1,12];
r:=rand(1..nops(A)):
seq(A[r()], i=1..3);

Since each time a number is chosen randomly from the entire list (a sample with a return), it can happen that the same number will be selected during the re-selection. It is also not difficult to implement a non-refund sample.

restart;
test := -2+exp(theta)+exp(-theta);
subs(exp(theta)=1/exp(-theta), test);
factor(%);

The final result:       (exp(-theta)-1)^2/exp(-theta)

Edit.

First write these data as a list of lists, then use the  plot  command as in this toy example:

XY:=[[1,2],[2,4],[3,5],[4,4]]:
plot(XY);

In the list of options for  pdsolve  command (with  numeric  option), the option  output = listprocedure  does not exist. See help on  ?pdsolve/numeric

Sol:=dsolve({diff(g(r),r,r)- r/R*g(r)=0, g(2*R)=0, D(g)(0)=R}) assuming R>0;
plot(eval(rhs(Sol),R=1), r=0..4, -2..2);

eq:=(-2*cos(x)^2+2*sin(x+(1/4)*Pi)^2-1)/sqrt(-x^2+4*x) = 0:
[solve(eq, x, allsolutions)];
subsindets(%, suffixed(_Z, integer), t->k);
map(t->unapply(t,k), %);
seq(`if`(evalf(p(0))>0 and evalf(p(0))<4,p(0),`if`(evalf(p(1))>0 and evalf(p(1))<4,p(1),NULL)), p=%);

  The final result:                      (1/2)*Pi,  (1/4)*Pi,  5*Pi*(1/4)

I remind that the original question was how to get a solution in the form  k*Pi and -Pi/48 + k*Pi/8
This form of the answer is based on a well-known identity  sin(a)-sin(b)=2*sin((a-b)/2)*cos((a+b)/2) . I don’t know a simple way to get this identity automatically in Maple, so I wrote it manually:

restart;
eq := sin(9*x-(1/3)*Pi) = sin(7*x-(1/3)*Pi):
a:=9*x-(1/3)*Pi:
b:=7*x-(1/3)*Pi:
solve(2*sin((a-b)/2)*cos((a+b)/2), allsolutions):
subsindets([%],suffixed(_Z, integer), t->k)[];

I replaced these assignments  x[0]:=1  and so on, due to which all your troubles:

Optimal_control_new.mw

plots:-implicitplot3d(`if`(y<x and z<y,cos(x-y)*cos(y-z)*cos(x-2*z)^3-0.6, NULL), x=0..5, y=0..5, z=0..2, style=surface, color=green, numpoints=500000, axes=normal);

remove~(is, S);

# Or (for old versions of Maple)

map(t->remove(is, t), S);