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

Example:

Points:=[[1,1],[1,3],[2,2],[3,3],[3,1]]:  # The points in the form of a list of their coordinates
P:=plot(Points, style=point, color=red, symbol=solidcircle, symbolsize=20):  # The points
L:=plot(Points, color=blue, thickness=2):  # The line segments
plots:-display(P, L, view=[0..4, 0..4]);  # All together

Edit.

with(plots):

K:=9;
deG:=diff(theta(t),t,t) + mu*diff(theta(t),t)+K*sin(theta(t))= 0;
deL:=diff(theta(t),t,t) + mu*diff(theta(t),t)+K*theta(t)= 0;
Iv:=theta(0)=0.75, D(theta)(0)=2.0;
dom1:=t=0..10;
soln1a:=dsolve({eval(deL,mu=0),Iv});

soln1b_1:=dsolve({eval(deL,mu=1),Iv});

gr1c:=multiple(plot,[rhs(soln1a),dom1,color=blue],[rhs(soln1b_1),dom1,color=purple]);
 

Addition.

In fact, there is no need to call plots package and use multiple command. It's simpler and smarter to use the standard  plot command, and write expressions to be plotted as a list:

restart;
K:=9;
deG:=diff(theta(t),t,t) + mu*diff(theta(t),t)+K*sin(theta(t))= 0;
deL:=diff(theta(t),t,t) + mu*diff(theta(t),t)+K*theta(t)= 0;
Iv:=theta(0)=0.75, D(theta)(0)=2.0;
dom1:=t=0..10;
soln1a:=dsolve({eval(deL,mu=0),Iv});
soln1b_1:=dsolve({eval(deL,mu=1),Iv});

plot([rhs(soln1a),rhs(soln1b_1)], dom1, color=[blue,purple]);  # Plotting by the usual  plot  command

Edit.