@PatrickT   I deleted the extra bracket.

Picture of output loaded in another way.

@Gauss   On my machine, the code is executed in a few seconds (in Classic M 12).  Run  restart  command  first.

@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)