@toandhsp  No, it can not. AreIsometric  procedure works only with the sets in the Euclidean plane, rather than in space.

@Alejandro Jakubi  Thanks for the explanation of the reason of the error.

@taro yamada  My method is only suitable for Classic Worksheet. If you write

A:=x/2+y/2:

normal(A/a);

then you see    

And then you just replace  symbol  a  by  the empty symbol  `` and get the desired output:

In the previous embodiment, the diagonal entries are always  0 , and should be random integers from  -9  to  9 . New version is free from this drawback:

roll:=rand(-9 .. 9):
Matrix(4, {seq((i, i)=roll(), i=1..4), seq(seq((i, j)=roll()+I*roll(), j=i+1..4), i=1..3)}, shape=hermitian);

Factoring  procedure takes 1/2 sec for this calculation, and the improvements made by Joe Riel and Carl Love even less time. See  http://www.mapleprimes.com/posts/141668-Partitions-Of-A-Natural-Number-Into-Factors

ts:=time():

Factoring(9!):

nops(%);

%%[1..20];

time()-ts;

I have already answered a similar question about the matrix  A . About the matrix  B  all is still easier. There is an inaccuracy in determining  B. Apparently instead of  -1/(2*M*(M-2))   should be -1/(M*(M-2)) , otherwise a contradiction with special cases.

@Bendesarts   Do all things similarly, only replace  eq  by  lhs(eq_liaison[1, 1])  and etc.

@Carl Love  Visually  f[`3x`,`4y`]  is identical to  f[`3x, 4y`]

@yihezhi 

f[`{6x}`]; f[`{3x,4y}`];

@Bendesarts  For some reason your last file is not downloaded. In my M 16 all works. See the attached file.

Solving_equation.mw

Here is the copy of the code:

eq := 2*sin(alpha(t))*sin(beta(t))*cos(a[1])*sin(a[1])*cos(gamma[1](t))*e[1]*rBTP[1]- 2*cos(alpha(t))*cos(beta(t))*h*h[1]+2*sin(beta(t))*sin(a[1])*h*rF[1]- 2*cos(alpha(t))*cos(a[1])^2*rBTP[1]*rF[1]-2*z(t)*sin(alpha(t))*cos(a[1])*rBTP[1]+ 2*z(t)*cos(alpha(t))*cos(beta(t))*h+2*sin(alpha(t))*cos(a[1])*h[1]*rBTP[1]+ 2*cos(beta(t))*cos(gamma[1](t))*e[1]*rBTP[1]+2*cos(beta(t))*cos(a[1])^2*rBTP[1]*rF[1]+ 2*sin(gamma[1](t))*e[1]*h[1]-2*z(t)*sin(gamma[1](t))*e[1]-2*cos(gamma[1](t))*e[1]*rF[1]- 2*cos(beta(t))*rBTP[1]*rF[1]-2*z(t)*h[1]-l[1]^2+ 2*sin(alpha(t))*cos(beta(t))*cos(a[1])*cos(gamma[1](t))*h*e[1]- 2*cos(alpha(t))*sin(beta(t))*sin(a[1])*sin(gamma[1](t))*e[1]*rBTP[1]- 2*sin(alpha(t))*sin(beta(t))*cos(a[1])*sin(a[1])*rBTP[1]*rF[1]+z(t)^2+ h[1]^2+rBTP[1]^2+e[1]^2+rF[1]^2+h^2+2*cos(alpha(t))*cos(a[1])^2* cos(gamma[1](t))*e[1]*rBTP[1]-2*cos(alpha(t))*sin(beta(t))*sin(a[1])*h[1]*rBTP[1]- 2*cos(alpha(t))*cos(beta(t))*sin(gamma[1](t))*h*e[1]-2*sin(beta(t))*sin(a[1])* cos(gamma[1](t))*h*e[1]-2*cos(beta(t))*cos(a[1])^2*cos(gamma[1](t))*e[1]*rBTP[1]+ 2*z(t)*cos(alpha(t))*sin(beta(t))*sin(a[1])*rBTP[1]-2*sin(alpha(t))*cos(beta(t))* cos(a[1])*h*rF[1]+2*sin(alpha(t))*cos(a[1])*sin(gamma[1](t))*e[1]*rBTP[1];

eq1 := collect(subs(sin(gamma[1](t)) = x, cos(gamma[1](t)) = y, eq), [x, y]);

A := coeff(eq1, x);

B := coeff(eq1, y);

C := tcoeff(eq1, [x, y]);

sol := solve(a*sin(t)+b*cos(t)+c, t);

gamma1 := unapply(eval(sol[1], {a = A, b = B, c = C}), t);

gamma2 := unapply(eval(sol[2], {a = A, b = B, c = C}), t);

diff(f(x,y),x$6);

diff(f(x,y),[x$3],[y$4]);

@ctnaz   No, the numbers  2  and  0.94  I got just by selection, banishing my program several times for different values​​.

@Preben Alsholm  Thanks for the solution! I've tried similar variants, but my mistake was that I took too little eps=10^(-6)

@Kitonum  

My first program was unsuccessful. I forgot that break keyword only interrupts the inner loop, but the program continues to run. Therefore returned not the first but the last solution.
One solution (the first one) can be found immediately when the program as a procedure to issue and uses  return keyword:

restart;
P := proc ()
local a, b, c, d, e, f, g, h;
for c from -10 to 10 do
for d from -10 to 10 do
for e from -10 to 10 do
for f from -10 to 10 do
for g from -10 to 10 do
for h from -10 to 10 do
a := c+2*d+4*e+6*f+10*g+22*h-16; b := -2*c-3*d-5*e-7*f-11*g-23*h+20;
if -10 <= a and a <= 10 and -10 <= b and b <= 10 then return [a, b, c, d, e, f, g, h]
end if end do end do end do end do end do end do
end proc:

P();

                                               [10, 10, -10, -10, -10, 8, 7, -1]

@Mikhail Drugov  Yes,  I think there is no automatic way.