You alone are meant to guarantee the applicability of the facility you are using. For linear algebra package, elements of the matrix should have a linear basis as arbitrary transformations required to obtain one will not be attempted. In your sample, the transformation sqrt(6)=sqrt(2)*sqrt(3) is not attempted and so the procedure is blind to existing 0s and it does not find the nullspace.

restart; with(linalg):
M:=matrix([[5/2, sqrt(3/2), sqrt(3/4)],
            [sqrt(3/2), 7/3, sqrt(1/18)],
            [sqrt(3/4), sqrt(1/18), 13/6]]):

eiv := proc(M)
    Normalizer := proc(x)
        local f, y;
        option `Shake things around.`;
            f := x -> normal(radnormal(x, rationalized)); y := radsimp(x, ratdenom); normal(f(x - y) + y)
        end;
    Testzero := proc(x) option `Shaken.`; evalb(Normalizer(x) = 0) end;
    [eigenvectors](map(Normalizer, M), radical)
end

E:=eiv(M);
P:=`@`(eval,augment,op,map2)(op@map,normalize,map2(op,3,E)):
map(radnormal,multiply(htranspose(P),M,P));

E:=[[4, 1, {}], [2, 1, {}], [1, 1, {}]]

Normalizer, Testzero are enviroment procedures (look up ?_Env). All modifications terminate with the scope.

You could recover a system by eliminating the RootOfs. The cheap way to do this is to use eliminate.

#interface(version);
#   Maple Worksheet Interface, Release 4, IBM INTEL NT, Apr 16 1996

recover := proc(S, T0, V0)
local Rs, n, Rs1, r, P, E, T, V;
option `This can fail in ways too many...`;
    Rs := indets(S, RootOf);
    Rs := remove(proc(x, S) has(op(1, x), S minus {x}) end, Rs, Rs);
    n := nops(Rs);
    if nargs < 2 then T := {} else T := T0 fi;
    if nargs < 3 then V := [] else V := V0 fi;
    if 0 < n then
        Rs1 := op(1, Rs);
        P := ({evala}@Norm)(r - Rs1);
        T := T union P;
        V := [op(V), Rs1 = r];
        E := eliminate(subs(op(-1, V), S) union P, r);
        if 1 < n then recover(E[2], T, V) else E[2], T, V fi
    else S, T, V
    fi
end;

Eliminate may return multiple results so you may want to fork it there.

_EnvExplicit:=false:
A:=evala({a = s/RootOf(_Z^2-s^2+s), b = -RootOf(_Z^2-s^2+s)/s, c = RootOf(_Z^2-s^2+s)});
rA:=recover(A);
rAs:=solve(rA[1],{a,b,c});
A2:=solve(rA[1],{a,b,s});
rA2:=recover(A2);
rA2s:=solve(rA2[1],{a,b,c});

A := {b = -RootOf(_Z^2-s^2+s)/s, c = RootOf(_Z^2-s^2+s), a = RootOf(_Z^2-s^2+s)/(s-1)}
rA := {c-a*s+a, s^2-s-a^2*s^2+2*a^2*s-a^2, b*s+a*s-a}, {-s^2+s+r^2}, [RootOf(_Z^2-s^2+s) = r]
rAs := {b = -RootOf((s-1)*_Z^2-s)*(s-1)/s, a = RootOf((s-1)*_Z^2-s), c = RootOf((s-1)*_Z^2-s)*(s-1)}
A2 := {a = RootOf(-c-_Z+_Z^2*c), s = c*RootOf(-c-_Z+_Z^2*c), b = -(c*RootOf(-c-_Z+_Z^2*c)-1)/c}
rA2 := {b*c+a*c-1, -c-a+a^2*c, -s+a*c}, {(-c-r+r^2*c)/c}, [RootOf(-c-_Z+_Z^2*c) = r]
rA2s := {b = -RootOf(_Z^2-s^2+s)/s, c = RootOf(_Z^2-s^2+s), a = RootOf(_Z^2-s^2+s)/(s-1)}

You appear to be solving for F(x,y) defined as y=F(x,y). So we can do the following:

EQ:=y=F(x,y);
sd:=map((f,x)->f=f(x),map2(`@@`,D,[$0..1])(y),x);
dEQ:=subs(sd,{EQ} union implicitdiff({EQ},{y}(x),{y},x));

EQ := y = F(x,y)
sd := [y = y(x), D(y) = D(y)(x)]
dEQ := {y(x) = F(x,y(x)), D(y)(x) = -D[1](F)(x,y(x))/(-1+D[2](F)(x,y(x)))}

So D(y)(x) = -D[1](F)(x,y(x))/(-1+D[2](F)(x,y(x))) contains the derivatives sought.

Using the definition F = min/max

minovermax:=piecewise(x<y,x,x=y,(x+y)/2,y)/piecewise(x<y,y,x=y,(x+y)/2,x):
rf:=map(`@`(unapply,L->(L[1],op(L[2])),[`@`(unapply(minovermax,x,y),op),L->L]),[assume(x0<y0,x1=y1,x2>y2),[x0,y0],[x1,y1],[x2,y2]]):
pf:=piecewise(x<y,rf[1](x,y),x=y,rf[2](x,y),x>y,rf[3](x,y)):
sol:=map2(proc(X,L,x,y) 'solve'({op(L[2],X)=y,op(L[1],X)},{x,y}) end,pf,[[1,2],[3,4],[5,6]],x,y);
sol;

sol := [solve({x < y, x/y = y},{x, y}), solve({1 = y, x = y},{x, y}), solve({y < x, y/x = y},{x, y})]
[{x = y^2, 0 < y, y < 1}, {y = 1, x = 1}, {y = 0, 0 < x}, {y < 1, x = 1}]

Above and to the right the x,y curve for {y=F(x,y), F=min/max}:

And a 3d display of y=F(x,y) in [x,y,F(x,y)] coordinates:

Plot1:=plots[display](
map(plot,[[y^2,y,y=0..1],[x,0,x=-2..0],[1,y,y=1..3]],color=grey,thickness=3),
plots[pointplot]([[1,1]],symbol=BOX,color=black),
axes=boxed
):
Plot2:=plots[display](
plottools[transform]((x,y)->[x,y,y])(Plot1),
labels=['x,y,f'],scaling=constrained):
plots[display](Plot2);

Since that curve contains your definitions then you should contain your differentiation to it.

Observe your system in detail first.

#Some Maple V code:

solveit := proc() frontend(solve, [args], [{Non(radical)}, {}]) end:

aaghulu := {-6-4*y-x-(1+y)*x+sqrt((4*(1+y))*(2+x)*(4+2*y+x)+(-(1+y)*x+2+x)^2), (2*(4+2*y+x))*(1+y)-(1+y)*x+2+x+sqrt((4*(1+y))*(2+x)*(4+2*y+x)+(-(1+y)*x+2+x)^2)-(2+y)*(-(1+y)*x+2+x+sqrt((4*(1+y))*(2+x)*(4+2*y+x)+(-(1+y)*x+2+x)^2))};usys:={6+2*x+4*y+y*x-u};
result:=`@`(factor,simplify)(aaghulu,usys,{x,y,u}) union usys;

tryit:=unapply('result union {(u^2)^(1/2)=someu},someu'):

solveit(tryit(u));
solveit(tryit(u),{sqrt(u^2),y});
solveit(tryit(-u),{sqrt(u^2),y,u});

#Results:

aaghulu := {2*(4+2*y+x)*(1+y)-(1+y)*x+2+x+((6+4*y+2*x+y*x)^2)^(1/2)-(2+y)*(-(1+y)*x+2+x+((6+4*y+2*x+y*x)^2)^(1/2)), -6-4*y-x-(1+y)*x+((6+4*y+2*x+y*x)^2)^(1/2)}
usys := {6+4*y+2*x+y*x-u}
result := {-u+(u^2)^(1/2), 6+4*y+2*x+y*x-u, (u-(u^2)^(1/2))*(1+y)}
{u = 6+4*y+2*x+y*x, x = x, (u^2)^(1/2) = 6+4*y+2*x+y*x, y = y}
{(u^2)^(1/2) = u, y = -(6+2*x-u)/(4+x)}
{u = 0, (u^2)^(1/2) = 0, y = -2*(3+x)/(4+x)}

For the solution of (u^2)^(1/2) = u if we assume u is real, positive and x,y are real then a map of x,y is:

plots[contourplot](6+4*y+2*x+y*x,x=-4..4,y=-4..4,contours=map(`*`,[$0..100],0.5));

Simply:

map(factor@simplify,tmp,map(numer@(lhs-rhs),{rule3}),{w[1],w[3],s[3]});

In general given polynomial reductions you should experiment with Maple's tools.

This function is probably sufficient for most expressions. If you are having malformed results then please post a response with an example.

function_coeffs := proc(A, v::set(name))
local S, T;
    S := indets(A, {function});
    S := select(has, S, v);
    T := {Non(map(identical, S))};
    frontend(proc(A, S) local V; [coeffs](A, S, 'V'), [V] end, [A, S union v], [T, {}])
end:

A:=a(x)*v*u+v*D(u)-D(v)*u:
function_coeffs(A,{u,v});

[-1, a(x), 1], [D(v)*u, u*v, v*D(u)]

Using eliminate we can get a result of

{S[1] = X*(Y*Q+sum(S[j],j = 2 .. i)*D*B+B*S[i+1]*D)/Y/(C*A+r[1]*X)}, {S[i+1] = RootOf(r[2]*_Z^2*X+(-A*S[i]*C+B*D*X+X*r[2]*Y)*_Z-A*S[i]*C*Y)}, {-A*S[i-1]*C*Y+A*S[i]*C*Y+B*S[i]*D*X+r[1]*S[i]*X*Y};

The third brace contains the remainder from eliminating S[1],S[i+1] from eq1,eq2,eq3.

You could try the following:

restart;
ohgamma:=gamma; #old constant gamma
alias(gamma=mygamma); #take over this presentation for your mygamma
{ohgamma,gamma,mygamma}; #multiplicity
solve(gamma^2+1,{gamma}); #now using mygamma guised as gamma
solve(gamma^2+1,{mygamma}); #now using mygamma and its guise
traperror(solve(gamma^2+1,{ohgamma})); #now using guise and old gamma

In your execution of the code did you get these values for s1,s2:

's1'=1/4*H*(exp(4*(4+C11)/H)-1)*exp(-2*(4+C11)/H)-1/4*H*(exp(4/H*C11)-1)*exp(-2/H*C11
's2'=-1/4*H*(exp(4*(20+C12)/H)-1)*exp(-2*(20+C12)/H)+1/4*H*(exp(4*(4+C12)/H)-1)*exp(-2*(4+C12)/H);

With these values I turned {r1,r2,r3,r4} into:

{cosh(1/H*C12) = RootOf(-H*exp(1/H)^8-2*C21*exp(1/H)^8+5*exp(1/H)^16-C22*exp(1/H)^16+C21*exp(1/H)^16-5+C21-C22+2*H*_Z^2)/exp(1/H)^4, sinh(1/H*C12) = 1/2*(20*exp(1/H)^48+4*exp(1/H)^48*C22-2*C21*exp(1/H)^8-2*C21*exp(1/H)^88+5*exp(1/H)^16+5*exp(1/H)^96-C22*exp(1/H)^16-C22*exp(1/H)^96+C21*exp(1/H)^16+C21*exp(1/H)^96-5-5*exp(1/H)^80+C21+C21*exp(1/H)^80-C22-C22*exp(1/H)^80)/RootOf(-H*exp(1/H)^8-2*C21*exp(1/H)^8+5*exp(1/H)^16-C22*exp(1/H)^16+C21*exp(1/H)^16-5+C21-C22+2*H*_Z^2)/exp(1/H)^4/H/(exp(1/H)^80-1), sinh(1/H*C11) = -(C21*exp(1/H)^80-exp(1/H)^72*C21+exp(1/H)^72*C22+5*exp(1/H)^72-10*exp(1/H)^40-2*C22*exp(1/H)^40-5*exp(1/H)^8-C21*exp(1/H)^8+C22*exp(1/H)^8+C21)/RootOf(2*H*_Z^2-H-2*C21)/H/(exp(1/H)^80-1), cosh(1/H*C11) = RootOf(2*H*_Z^2-H-2*C21)};

Final solution used RootOfs for degree 4 polynomials.

The following code shows the two equations, here in variables A, B and constants alpha, beta, are linear transforms of same as captured in: h(alpha*A,beta) = g(alpha*A,B), h(A,beta) = g(A,B). Relations are stored in EQ1s, EQf and EQ1.

restart;
eq1:=(1-Pi*x/2)*(f/2/E)*yc-1/3=Int(t^4/sqrt((t^2-(f/2/E)*yc)^2+ya^2),t=0..1);
eq2:=x^2*yc*(1-Pi*x/2)-1/3=Int(t^4/sqrt((t^2-x^2*yc)^2+ya^2),t=0..1);

eq2:=subs(t^4-2*t^2*x^2*yc+x^4*yc^2+ya^2 = (t^2-x^2*yc)^2+ya^2,eq2);

EQ1s:=A = 1/2*yc/E*f, B = ya^2, alpha = 2*E*x^2/f, beta = x;
G:=(q,r)->Int(t^4/((t^2-q)^2+r)^(1/2),t = 0 .. 1);;
H:=(q,r)->q*(1-1/2*Pi*r)-1/3;
EQf:=h=H,g=G;
EQ1:=h(alpha*A,beta) = g(alpha*A,B), h(A,beta) = g(A,B);

#show
g=eval(G);h=eval(H);
EQf;
EQ1s;
EQ1;
form:=(eval@subs)([EQf],[EQ1s],[EQ1]);
test:=map(evalb,form-[eq2,eq1]);
(eval@subs)(EQf,[EQ1]);
EQ1;

The matrix M can be inverted through the use of diagonalization. The inverse matrix entries are very large. Direct computation of the inverse fails with object too large error. A proposition of computation path follows. Observed sturcture of M matrix is posted in M_Qver (using variables Q1..Q7) and observed sturcture of M inverse matrix is posted in Mi_Pver (using variables P1..P12). You can solve for the new variables.
Use equation
M = M_Qver  
to get Q variables from a,b,c,d,e,w. Convert
0 = M_Qver &* Mi_Pver - 1.
0 = Mi_Pver &* M_Qver - 1
into a system of linears in variables P1..P12. to get P variables from Q variables. Solve for P variables and substitute into Mi_Pver the inverse of M. In following code matrix PMi last column contains values of P1..12. This computes in minutes.

Code:

restart:
with(linalg):

M:=matrix([ [a+b+c+d+e,a+b+d+w*e,a+b+d+w^2*e,a+b+c,a+b,a+b,a+c,a, a,a+c,a,a], [a+b+d+w*e,a+b+c+d+e,a+b+d+w*e,a+b,a+b+c,a+b,a,a+c,a, a,a+c, a], [a+b+d+w^2*e,a+b+d+w*e,a+b+c+d+e,a+b,a+b,a+b+c,a,a,a+c,a, a,a+c], [a+b+c,a+b,a+b,a+b+c+d+e,a+b+d+w*e,a+b+d+w^2*e,a+c,a,a, a+c,a, a], [a+b,a+b+c,a+b,a+b+d+w*e,a+b+c+d+e,a+b+d+w*e,a,a+c,a,a, a+c, a], [a+b,a+b,a+b+c, a+b+d+w^2*e,a+b+d+w*e, a+b+c+d+e,a,a,a+c, a, a,a+c], [a+c, a,a,a+c,a, a,a+b+c+d+e,a+b+d+w*e,a+b+d+w^2*e, a+b+c, a+b,a+b], [a,a+c,a,a, a+c,a,a+b+d+w*e,a+b+c+d+e,a+b+d+w*e, a+b, a+b+c, a+b], [a,a, a+c,a, a, a+c, a+b+d+w^2*e,a+b+d+w*e,a+b+c+d+e, a+b,a+b,a+b+c], [a+c,a,a,a+c,a,a,a+b+c, a+b,a+b,a+b+c+d+e, a+b+d+w*e, a+b+d+w^2*e], [a, a+c,a,a,a+c,a,a+b,a+b+c, a+b, a+b+d+w*e, a+b+c+d+e, a+b+d+w*e], [a,a,a+c,a,a,a+c,a+b, a+b,a+b+c,a+b+d+w^2*e,a+b+d+w*e, a+b+c+d+e]]);

M_subs:={a = Q3, a+b = Q6, a+b+c = Q7, a+b+d+w^2*e = Q2, a+b+d+w*e = Q5, a+b+c+d+e = Q1, a+c = Q4};
Q_vals := {Q1 = a+b+c+d+e, Q5 = a+b+d+w*e, Q7 = a+b+c, Q2 = a+b+d+w^2*e, Q4 = a+c, Q3 = a, Q6 = a+b};
Q_ident:={Q3 = Q4-Q7+Q6};

M_Qver := matrix([[Q1, Q5, Q2, Q7, Q6, Q6, Q4, Q3, Q3, Q4, Q3, Q3], [Q5, Q1, Q5, Q6, Q7, Q6, Q3, Q4, Q3, Q3, Q4, Q3], [Q2, Q5, Q1, Q6, Q6, Q7, Q3, Q3, Q4, Q3, Q3, Q4], [Q7, Q6, Q6, Q1, Q5, Q2, Q4, Q3, Q3, Q4, Q3, Q3], [Q6, Q7, Q6, Q5, Q1, Q5, Q3, Q4, Q3, Q3, Q4, Q3], [Q6, Q6, Q7, Q2, Q5, Q1, Q3, Q3, Q4, Q3, Q3, Q4], [Q4, Q3, Q3, Q4, Q3, Q3, Q1, Q5, Q2, Q7, Q6, Q6], [Q3, Q4, Q3, Q3, Q4, Q3, Q5, Q1, Q5, Q6, Q7, Q6], [Q3, Q3, Q4, Q3, Q3, Q4, Q2, Q5, Q1, Q6, Q6, Q7], [Q4, Q3, Q3, Q4, Q3, Q3, Q7, Q6, Q6, Q1, Q5, Q2], [Q3, Q4, Q3, Q3, Q4, Q3, Q6, Q7, Q6, Q5, Q1, Q5], [Q3, Q3, Q4, Q3, Q3, Q4, Q6, Q6, Q7, Q2, Q5, Q1]]);

Mi_Pver := matrix([[P10, P9, P1, P2, P5, P3, P6, P7, P8, P6, P7, P8], [P9, P11, P9, P5, P12, P5, P7, P4, P7, P7, P4, P7], [P1, P9, P10, P3, P5, P2, P8, P7, P6, P8, P7, P6], [P2, P5, P3, P10, P9, P1, P6, P7, P8, P6, P7, P8], [P5, P12, P5, P9, P11, P9, P7, P4, P7, P7, P4, P7], [P3, P5, P2, P1, P9, P10, P8, P7, P6, P8, P7, P6], [P6, P7, P8, P6, P7, P8, P10, P9, P1, P2, P5, P3], [P7, P4, P7, P7, P4, P7, P9, P11, P9, P5, P12, P5], [P8, P7, P6, P8, P7, P6, P1, P9, P10, P3, P5, P2], [P6, P7, P8, P6, P7, P8, P2, P5, P3, P10, P9, P1], [P7, P4, P7, P7, P4, P7, P5, P12, P5, P9, P11, P9], [P8, P7, P6, P8, P7, P6, P3, P5, P2, P1, P9, P10]]);

P_vector:=[P1, P2, P3, P4, P5, P6, P7, P8, P9, P10, P11, P12];
QP_set:=convert(evalm(M_Qver &* Mi_Pver - 1),set) union convert(evalm(Mi_Pver &* M_Qver - 1),set);
QP_set:=subs(Q_ident,QP_set):

PM:=genmatrix(QP_set,P_vector,true);
PMi:=gaussjord(PM,coldim(PM)-1):

Expecting real values over the real range you'd choose this function:

piecewise(0 < x,2*(-1+x)/x^(5/3),x = 0,FAIL,x < 0,(-1+I*3^(1/2))*(-1+x)/x^(5/3))

note this function is real valued but floating point valuation produces an imaginary component, so use Re.

plot(Re(piecewise(0 < x,2*(-1+x)/x^(5/3),x = 0,FAIL,x < 0,(-1+I*3^(1/2))*(-1+x)/x^(5/3))),x=-4..4,view=[-4..4,-4..4]);

or

plot(piecewise(0 < x,2*(-1+x)/x^(5/3),x = 0,FAIL,x < 0,-2*(-1+x)/(-x)^(5/3)),x=-4..4,view=[-4..4,-4..4]);

This was written with an old version of Maple so I apologize for any difficulties.

restart;
h:=x->-2/x^(5/3)+2/x^(2/3);
assume(x,complex,q,complex);
h(x);
rationalize(h(x));
E:=(denom@rationalize)(1/(q-rationalize(h(x))));
sol:=[solve](E,q);

atry:=map(x->piecewise(Im(x)=0,Re(x),FAIL),sol):
ans:=unapply(eval(atry,1),x)(-1/2);
ans:=remove(type,ans,map(identical,{FAIL}));
ans[1];
evalf(ans[1],20);

                                   2/3
                                6 2


                        9.5244063118091968488