Here is a slightly different way (more programmatic), which does not require you to manually determine the required positions in res :

restart;
 res := dsolve({25*(diff(y(t), t, t))+4*(diff(y(t), t))-3*y(t) = cos(3*t), y(0) = 0, D(y)(0) = 1}, numeric):
with(Optimization):
Maximize(s->eval(y(t), res(s)), 0..4);
Maximize(s->eval(diff(y(t),t), res(s)), 0..4);
 

Here is another way, based on the direct construction of a spherical cap. With this method, the cap is shown as a part of a sphere with the same gridlines on the surface. The names of the parameters in the procedure and their meaning are the same as in Rouben's one:

restart;
SphericalCap:=proc(cap_extent, colatitude, longitude)
local alpha, phi0, theta0, V0, V, L0, L, Var, A, B, C, T;
uses plots;
alpha:=cap_extent;  phi0:=longitude;  theta0:=colatitude;
V0:=<cos(phi0)*sin(theta0), sin(phi0)*sin(theta0), cos(theta0)>;
V:=<cos(phi)*sin(theta), sin(phi)*sin(theta), cos(theta)>;
L0:=convert(V0, list); L:=convert(V, list); Var:=<x,y,z>;
A:=plot3d(map(p->`if`(arccos(V0.V)<=alpha, p, NULL), L), phi=0..2*Pi, theta=0..Pi, color=yellow, numpoints=90000):
B:=plot3d(L, phi=0..2*Pi, theta=0..Pi, numpoints=90000):
C:=plots:-intersectplot(surface(V0.(Var-cos(alpha)*V0), x=-1..1, y=-1..1, z=-1..1), surface(L, phi=0..2*Pi, theta=0..Pi), color=black, thickness=2):
T:=plots:-textplot3d([[1.37,0,0, x], [0,1.37,0, y], [0,0,1.37, z]], align=[right,below], font=[times,bold,18]):
plots:-display(A,B,C,T, axes=normal, view=[-1.37..1.37,-1.37..1.37,-1.37..1.37], orientation=[50,65]);
end proc:

Example of use:

SphericalCap(Pi/6, Pi/6, Pi/2);

SphericalCap.mw

For example, here is the solution for A:
solA:=solve(Determinant(Matrix([[xx1,yy1,1],[xx2,yy2,1],[xx3,yy3,1]])) =(1/2)*aa*d*s*u+(1/2)*aa*d*s*a*t+(1/2)*d*v*u*t+(1/4)*d*v*a*t^2, {xx1,yy1,xx2,yy2,xx3,yy3});

As the values for  yy1, xx2, yy2, xx3, yy3  , you can take any numbers, only should be  yy2<>yy3

For example:
map(t->lhs(t)=eval(rhs(t), [xx2 = 1, xx3 = 2, yy1 = 3, yy2 = 5, yy3 = 4]), solA);

Edit.

Your code is incomplete. I added the initial conditions (arbitrarily) and something else:

restart;
M:=<1,0; 0,1>: d2y:=<diff(y__1(t),t,t),diff(y__2(t),t,t)>: K:=<3,1; 1,6>:
y:=<y__1(t),y__2(t)>: F:=<10,5>:
eq:=M.d2y+K.y=F:
ic:=y__1(0)=1, y__2(0)=2, D(y__1)(0)=0, D(y__2)(0)=-1:
A:=(lhs-rhs)(eq);
sol:=dsolve({convert(A,list)[ ], ic}, numeric);
plots:-odeplot(sol, [[t,y__1(t)],[t,y__2(t)]], t=0..10, color=[red,green]);
solve(A[1], diff(y__1(t),t,t));
plots:-odeplot(sol, [t,%], t=0..10, color=blue, labels=[t,diff(y__1(t),t,t)]);
 

Addition: hereafter send your code in text form (1d math) so that everyone can simply copy it to a worksheet. Therefore, I recommend always type the code in 1d math (as I did above). This is faster and allows you to better understand the Maple syntax.

Edit.

Add the following line to your code:
map(simplify~, %[2], zero);

See corrected file:

ANALYTIC_1_1.mw

Do this in polar coordinates:

plot3d(eval([x, y, x^3-2*x^2*y], [x=r*cos(phi), y=r*sin(phi)]), r=0..1, phi=0..2*Pi, axes=normal, numpoints=5000);

In Cartesian coordinates, this can also be done, but the quality of plotting is slightly worse:

plot3d(x^3-2*x^2*y, x=-1..1, y=-sqrt(1-x^2)..sqrt(1-x^2), axes=normal, numpoints=5000);

Edit.
 

You have come up with an intricate and very inefficient method of solving the problem, in which you made a lot of mistakes. Such tasks are instantly solved using seq command:

make_PSI:=nplanets->[seq(2*Pi*k/nplanets, k=0..nplanets-1)]:

Example of use:

make_PSI(4);

                             [0, (1/2)*Pi, Pi, 3*Pi*(1/2)]
 

with(plots):
Sys:=[x(t) = cos(t) + t*sin(t), y(t) = sin(t) - t*cos(t)]:
C:=spacecurve([rhs~(Sys)[ ], 0], t=0..Pi/2, axes=normal, color=red, thickness=4): # The curve in 3d
S:=plot3d(eval([x(t), y(t)*cos(phi), y(t)*sin(phi)], Sys), phi=0..2*Pi, t=0..Pi/2): # The surface of rotation
display(C, S, style=surface, scaling=constrained, labels=[x,y,z], orientation=[70,70]);  # Alltogether

Addition - animation of the rotation:

RC:=phi->plottools:-rotate(C, phi, [[0,0,0],[1,0,0]]):
RS:=phi->plot3d(eval([x(t),y(t)*cos(s),y(t)*sin(s)],Sys), s=0..phi, t=0..Pi/2, style=surface):
animate(display, ['RC'(phi), 'RS'(phi)], phi=0..2*Pi, scaling=constrained, orientation=[70,70], lightmodel=light1, frames=90);

SurfaceOfRotation.mw

eq := diff(x(t), t) = x(t) + 1:
indets(eq, {name, function(name)});

                                {t, x(t)}

Maybe this is a bug (the first version with the loop)?

Here is a workaround:

nonIdMaps := [ ]:
F:=[x -> x,x -> 2*x,x -> 3*x]:
for i from 1 to 3 do
F[i], F[i](y);
if F[i](y) <> y then nonIdMaps := [nonIdMaps[],F[i]] end if;
end do;
nonIdMaps, map(f -> f(y),nonIdMaps);

Try

labels = [x, `#mover(mi("sigma"),mo("&circ;"))`[y]]
 

# Or

labels = [x, conjugate(sigma[y])]

Edit.

See help on  mtaylor  command.

Example. Should be an initial condition for the numerical solution:

sol:=dsolve({diff(y(x),x) = -2*x*y(x) + 1, y(0)=2}, numeric, method=classical, stepsize=0.1); # The numerical solution in the form of a procedure 
plots:-odeplot(sol, [x,y(x)], x=0..1); # The plotting of the solution in the specific range 

# Finding the values of the function y(x) at individual points
eval(y(x), sol(0.1));
eval(y(x), sol(0.15));
eval(y(x), sol(0.2));

We see that to find the values of the function at intermediate points Maple uses linear interpolation in this method.

dsolve.mw

Edit.

See  desired_form_in_output_new.mw