Or

sol:=dsolve(diff(y(x),x)= x/(sqrt(x^2-16))*1/(2*y(x)), implicit);
applyop(simplify, [1, 2], %);
solve(%, y(x));

or

sol:=dsolve(diff(y(x),x)= x/(sqrt(x^2-16))*1/(2*y(x)), implicit);
evalindets(%, `*`, simplify);
solve(%, y(x));

or

sol:=dsolve(diff(y(x),x)= x/(sqrt(x^2-16))*1/(2*y(x)), implicit);
solve(map(simplify, lhs(%)), y(x));
 

The procedure classifies the terms of an arbitrary multivariate polynomial with numeric coefficients according to the indicated characteristics:

ClassTerms:=proc(P)
local S, C, t, T, Terms, Tp, Tn, Tind;
S:=indets(P);
C:=[coeffs(P,S,'t')];
T:=[t];
Terms:=convert(zip(`[]`,C,T),set);
Tp:=select(p->p[1]>0 and `and`(seq(type(degree(p[2],s),even), s=S)), Terms);
Tn:=select(p->p[1]<0 and `and`(seq(type(degree(p[2],s),even), s=S)), Terms minus Tp);
Tind:=Terms minus `union`(Tp,Tn);
print(P_positive={seq(p[1]*p[2],p=Tp)});
print(P_negative={seq(p[1]*p[2],p=Tn)});
print(P_indeterminate={seq(p[1]*p[2],p=Tind)});
end proc:


Example of use:

P:=randpoly([x,y,z], dense, degree = 4);

ClassTerms(P);

ClassTerms.mw

If the function is linear, then you can delete this option in your code or do not use it at all if you make a linearity check as follows:

type(a*x+b*y+c, linear(x, y));
                              true 

type(a*x^2+b*y+c, linear(x, y));
                              false

Do

labels=[z,`&Delta;z`/`&Delta;t`]
 

See the toy example:

f:=x[1]->x[1]^2;

whattype(x[1]);

                       indexed

A workaround:

f:=x1 -> x1^2;  # Or
f:=x__1 -> x__1^2;

Two triangles  ABC  and  A'B'C'  are similar. The first triangle can be obtained from the second one by successive transformations: rotation, translation and stretching. Animation of this:
restart;
with(plots): with(plottools):
T1:=display(polygon([[-6,0],[32/3,0],[0,8]],color=cyan), textplot([[-6,0,A],[32/3,0,C],[0,8,B]], font=[times,bold,18], align=right), scaling=constrained, axes=none):
T2:=display(polygon([[14,0],[17,0],[17,-4]],color=yellow), textplot([[14,0,"A'"],[17,0,"B'"],[17,-4,"C'"]], font=[times,bold,18], align=right), scaling=constrained, axes=none):
f:=(p,phi)->rotate(p,phi,[14,0]):
An1:=animate(display,[f(T2,phi)],phi=0..arctan(4/3), background=T1):
T3:=f(T2,arctan(4/3)):
g:=(p,h)->translate(p,h,0):
An2:=animate(display,[g(T3,h)],h=0..-15.8, background=T1):
T4:=g(T3,-15.8,0):
h:=(p,k)->homothety(p,k):
An3:=animate(display,['h'(T4,k)],k=1..10/3, background=T1):
display([An1,An2,An3], insequence=true, size=[800,800]);

Example:

sol:=RootFinding:-Isolate(x^4-3*x^2+2);
eval(x, sol[1]);  # The first root

In your for loop  m  should be a numeric not a symbol.
In addition, if you want to refer to the entries of your matrix 2 by 2 in a for loop, then use a nested loop:
...
for i from 1 to 2 do
for j from 1 to 2 do
...

To find positive solutions to this system with some additional conditions, you do not need to solve it several times in a loop. Maple actually everything has already solved  for your assumption  assume(0<d, d<1). You must only substitute the desired values of the parameter _Z1 into the found formula for d. See

restart;
assume(0<d, d<1):
assume(-0.01<a, a<0):
sys:={Re((-80*Pi*I*a/((a+1)^3))*exp(4*Pi*I*d)) = -0.4, Im((-80*Pi*I*a/((a+1)^3))*exp(4*Pi*I*d)) = 0.8}:
sol:=solve(sys, {a,d}, useassumptions = true,AllSolutions=true);
about(_Z1);
seq(eval(sol[2],_Z1=z), z in [0,1]);

We found the first and second positive values  d  satisfying conditions  d>0  and  d<1

Maple correctly decides your system. See
restart;
assume(d::real, d>0):
assume(a::real, -0.01 < a, a < 0):
sys:={-800*Pi*a*cos(6.557377048*Pi*(3.470797713+d))/(a+1)^3 = -.9396060697, 800*Pi*a*sin(6.557377048*Pi*(3.470797713+d))/(a+1)^3 = -.3238482794};
solve(sys, {a,d}, useassumptions=true, AllSolutions=true);
about(_Z1);

We see that the system has an infinite number of solutions that satisfy the given conditions, and the smallest positive  d  will be for  _Z1=11

Visualization:

sol:=solve(sys, {a,d}, useassumptions=true, AllSolutions=true):
eval(rhs(sol[2]), _Z1=11);  # The smallest positive  d
plots:-implicitplot(convert(sys, list), a=-0.01..0, d=0..0.03, color=[red,blue], gridrefine=5); 

Edit.

restart;
convert(D(s)(t), diff);

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.