M. Hirnyk wrote "It is clear that the set of the real solutions of the equation under consideration is infinite". This is true, but it does not follow from the plot, since any plot is constructed from a finite number of points. For example, try using the plots to determine how many real roots have the following two equations:

1)  10^(10*x)-x^(10^10)+1 = 0

2)  log[0.0659878](x) = 0.0659878^x

For the solution of the original equation, we note that the critical points of the function  

x->cos(x)^2+2*x*cos(x)-2*sin(x)+4  

going by equal distances to each other excepting the unique negative root of the equation  x+cos(x)=0 :

restart;

A:=cos(x)^2+2*x*cos(x)-2*sin(x)+4:

B:=factor(diff(cos(x)^2+2*x*cos(x)-2*sin(x)+4, x));

_EnvAllSolutions := true: solve(B);

simplify(eval(A, x=Pi*k)) assuming k::posint;

We see that the signs of the function at these critical points are alternated, and because between these points the function is monotonic, then in each interval  x = Pi*(n-1) .. Pi*n ,  n>=1,  there is only one real root. For negative  x the situation is similar.

The procedure  nthRoot  finds the  n_th  real positive or negative root with  N  significant digits:

nthRoot:=proc(n, N, t)

Digits:=N;

if t=positive then fsolve(cos(x)^2+2*x*cos(x)-2*sin(x)+4, x=Pi*(n-1)..Pi*n)  else

fsolve(cos(x)^2+2*x*cos(x)-2*sin(x)+4, x=-Pi*(n+1)..-Pi*n)  fi;

end:

Example (the first 10 positive and negative roots):

seq(nthRoot(n, 15, positive), n=1..10);

seq(nthRoot(n, 15, negative), n=1..10);

eq1:=diff(u[1](z),z)=a[11](z)*u[1](z)+a[12](z)*u[2](z):

u[1]:=z->sin(z)*f(z):

eq1;

or 

eq1:=diff(u[1](z),z)=a[11](z)*u[1](z)+a[12](z)*u[2](z):

eval(eq1, u[1]=sin(z)*f(z));

Just one region occludes the other.

Workaround:

with(plots):

I1 := inequal({x < 1, y < 1}, x = 0 .. 1, y = 0 .. 1, optionsfeasible = [color = "Brown", thickness = 2], view = [0 .. 2, 0 .. 2]):

I2 := inequal({y > 1, x < 1}, x = 0 .. 1, y = 1 .. 2, optionsfeasible = [color = "Yellow", thickness = 2], view = [0 .. 2, 0 .. 2]):

display(I1);

display(I2);

display({I1, I2});

restart;

A:=(1/a[x])^(-k-1+epsilon)/(1/a[I])^(-k-1+epsilon);

Rule:=(x::anything)^(z::anything)/(y::anything)^(z::anything)=(x/y)^z:

B:=applyrule(Rule, A);

subs({op(1,B)=op(1,B)^(-1), op(2,B)=-op(2,B)}, B);  # Result

For example

simplify((a+b)^2, {a^2+b^2=1, a*b=1/2});

                                 2

Addition. The above example is trivial. Here is less obvious example:

simplify(a^20+b^20, {a+b=5, a*b=1});

                       40648568638127

Idea of ​​the proof:
If  x*(x +1)  is divisible by 100!  then because the numbers  x  and  x+1  are relatively prime,  in the expansion of  100!  into prime factors any  p^r  can belong to only one of the factors  x  or  x+1.

restart;

A:={op(ifactors(100!)[2])}:

B:=[seq(op(combinat[choose](A,k)), k=1..7)]:

C:=[seq([B[i], A minus B[i]], i=1..nops(B))]:

for c in C do

a:=mul(c[1,i,1]^c[1,i,2], i=1..nops(c[1])): b:=mul(c[2,i,1]^c[2,i,2], i=1..nops(c[2])):

isolve(a*u-b*v=1): x:=tcoeff(b*rhs(%[2])): if x<4*99! then break fi:

isolve(b*v-a*u=1): x:=tcoeff(a*rhs(%[1])): if x<4*99! then break fi:

od:

x;

Verification:

is(irem(x*(x+1), 100!)=0 and x<4*99!);

                               true

This method lets you easily find the smallest  x  satisfying the given conditions, but it will take much more time.

You can use  plot3d  command:

A:=mtaylor(sin(x*y), [x=1,y=2], 6);

plot3d(A, x=0..3, y=0..3, view=-2..2, axes=normal, scaling=constrained);

Here is an independent solution of the logical problem. Code is quite complicated, because based on the phase-solving, but it works fast. Carl, do  you get the same answers in your package?

restart;

Position:=[$ 1..5]:

Names:= [Bob, Keeley, Rachael, Eilish, Amy]:

TV:=[Simpsons, Coronation, Eastenders, Desperate, Neighbours]:

Dest:= [Fra, Aus, Eng, Afr,Ita]:

Ages:= [14, 21, 46, 52, 81]:

Hair:=[afro, long, straight, curly , bald]:

Lives:= [town, city, village, farm, youth]:

P:=proc(A,B)

if (A implies B) and (B implies A) then true else false fi;

end:

L:=[]:

for i in Position do

for j in Names do

for k in TV do

for p in Dest do

for q in Ages do

for s in Hair do

for t in Lives do

if convert([P(i=3,k=Desperate), P(i=1,j=Bob), (k=Simpsons implies t<>youth) and (t=youth implies k<>Simpsons), (p=Afr implies j<>Rachael) and (j=Rachael implies p<>Afr), P(t=village,q=52), P(p=Aus,s=straight), P(p=Afr,k=Desperate), P(q=14,i=5), P(j=Amy, k=Eastenders), P(p=Ita,s=long), P(j=Keeley,t=village), P(q=46,s=bald), P(i=4, p=Eng),  (k=Coronation implies s<>afro) and (s=afro implies k<>Coronation), (j=Rachael implies s<>afro) and (s=afro implies j<>Rachael), P(q=21,t=youth), P(k=Coronation, s=long), P(q=81,t=farm), P(p=Fra,t=town), (j=Eilish implies s<>straight) and (s=straight implies j<>Eilish)], `and`) then L:=[op(L), [i,j,k,p,q,s,t]] fi:

od: od: od: od: od: od: od:

M:=[ListTools[Categorize]((x,y)->x[1]=y[1], L)]:

R:=[]:

for i in M[1] do

for j in M[2] do

for k in M[3] do

for p in M[4] do

for q in M[5] do

if ((i[3]=Simpsons and j[7]=youth) or (j[3]=Simpsons and (i[7]=youth or k[7]=youth)) or (k[3]=Simpsons and (j[7]=youth or p[7]=youth)) or (p[3]=Simpsons and (k[7]=youth or q[7]=youth)) or (q[3]=Simpsons and p[7]=youth)) and ((i[2]=Rachael and (j[4]=Afr or k[4]=Afr or p[4]=Afr or q[4]=Afr)) or (j[2]=Rachael and (k[4]=Afr or p[4]=Afr or q[4]=Afr)) or (k[2]=Rachael and (p[4]=Afr or q[4]=Afr)) or (p[2]=Rachael and q[4]=Afr)) and ((i[3]=Desperate and j[3]=Neighbours) or (j[3]=Desperate and (i[3]=Neighbours or k[3]=Neighbours)) or (k[3]=Desperate and (j[3]=Neighbours or p[3]=Neighbours)) or (p[3]=Desperate and (k[3]=Neighbours or q[3]=Neighbours)) or (q[3]=Desperate and p[3]=Neighbours)) and ((i[3]=Coronation and j[6]=afro) or (j[3]=Coronation and (i[6]=afro or k[6]=afro)) or (k[3]=Coronation and (j[6]=afro or p[6]=afro)) or (p[3]=Coronation and (k[6]=afro or q[6]=afro)) or (q[3]=Coronation and p[6]=afro)) and ((i[2]=Rachael and j[6]=afro) or (j[2]=Rachael and (i[6]=afro or k[6]=afro)) or (k[2]=Rachael and (j[6]=afro or p[6]=afro)) or (p[2]=Rachael and (k[6]=afro or q[6]=afro)) or (q[2]=Rachael and p[6]=afro)) and not ((i[2]=Eilish and j[6]=straight) or (j[2]=Eilish and (i[6]=straight or k[6]=straight)) or (k[2]=Eilish and (j[6]=straight or p[6]=straight)) or (p[2]=Eilish and (k[6]=straight or q[6]=straight)) or (q[2]=Eilish and p[6]=straight))  then

R:=[op(R), [i,j,k,p,q]] fi:

od: od: od: od: od:

Result:=select(x->convert([seq(nops({seq(x[i,j], i=1..5)})=5, j=1..7)],`and`), R):

print(`First solution`);

Matrix(Result[1])^+;

print(``);

print(`Second solution`);

Matrix(Result[2])^+; 

Puzzle.mws

with(LinearAlgebra):

a, b, c:=<22,11,5>, <13,6,3>, <-5,-2,-1>:

a1:=a/norm(a,2):

b0:=b-(a1.b)*a1:

b1:=b0/norm(b0,2):

c0:=c-((a1.c)*a1+(b1.c)*b1):

c1:=c0/norm(c0,2):

a1,  b1,  c1;

Vectors  a1, b1, c1  are unit and mutually perpendicular:

a1.a1, b1.b1, c1.c1, a1.b1, a1.c1, b1.c1;

            1, 1, 1, 0, 0, 0

convert(op(0.25)[1], symbol)/(10^(-op(0.25)[2]));

 ``(op(0.25)[1])/``(10^(-op(0.25)[2]));

The second method allows you to use this fraction in subsequent calculations:

expand(%*3/5);

              3/20

Solution without  dsolve (only for rectangular tables):

restart;

A:=plots[animate](plot,[[2/Pi*abs(arcsin(sin(Pi/2*(0.5+0.2*t)))), 2/Pi*abs(arcsin(sin(Pi/2*(0.5+sqrt(3)/5*t)))), t=0..a], thickness=5, color=red, numpoints=1000], a=0..20, frames=200):

B:=plottools[rectangle]([0,1], [1,0], thickness=2, color=green):

plots[display](A, B, scaling=constrained, axes=none);

with(RealDomain):

 plot(x^(1/3), x=-8 .. 8, thickness=5, scaling=constrained);

For Classic M 12:

A:=1/2+(1/3)*sqrt(2);

numer(A)/(``*denom(A));

For Standard:

A := 1/2+(1/3)*sqrt(2);

numer(A)/convert(denom(A), symbol);

Another variant:

A := 1/2+(1/3)*sqrt(2);

``(numer(A))/denom(A);

For your example it's more reliable to use a linear model, which is obtained if to apply  ln  function to your equation:

ln(x)=ln(a)+b*ln(y)+c*ln(z)  or  x1=a1+b*y1+c*z1 , where  x1=ln(x) ,  a1=ln(a) ,  y1=ln(y) ,  z1=ln(z) 

xlist := [0.1e-1, 0.1e-1, 0.1e-1, 0.113e-2, 0.144e-2, 0.171e-2, 0.21e-2, 0.23e-2, 0.183e-2, 0.233e-2, 0.28e-2, 0.331e-2, 0.409e-2, 0.304e-2, 0.377e-2, 0.459e-2, 0.576e-2, 0.708e-2, 0.424e-2, 0.502e-2, 0.623e-2, 0.763e-2, 0.918e-2, 0.529e-2, 0.63e-2, 0.778e-2, 0.965e-2, 0.116e-1, 0.63e-2, 0.747e-2, 0.899e-2, 0.113e-1, 0.137e-1, 0.724e-2, 0.872e-2, 0.105e-1, 0.13e-1, 0.158e-1, 0.848e-2, 0.1e-1, 0.122e-1, 0.15e-1, 0.182e-1]:

ylist := [16.529, 13.333, 10.929, 2.66, 2.75, 2.77, 2.69, 2.58, 4.4, 4.67, 4.62, 4.51, 4.48, 7.5, 7.58, 7.64, 7.75, 7.77, 10.1, 10.1, 10.1, 10.2, 10.2, 12.6, 12.7, 12.8, 12.7, 12.8, 14.9, 15, 15, 15, 15, 17.3, 17.5, 17.4, 17.3, 17.3, 20, 19.9, 19.9, 19.9, 20]:

zlist := [100, 1000, 10000, 1, 10, 100, 1000, 10000, 1, 10, 100, 1000, 10000, 1, 10, 100, 1000, 10000, 1, 10, 100, 1000, 10000, 1, 10, 100, 1000, 10000, 1, 10, 100, 1000, 10000, 1, 10, 100, 1000, 10000, 1, 10, 100, 1000, 10000]: 

X:=convert(map(ln,xlist), Vector):

Y:=convert(map(ln,ylist), Vector):

Z:=convert(map(ln,zlist), Vector):

YZ:=Matrix([Y,Z]):

Statistics[LinearFit]([1,y1,z1], YZ, X, [y1,z1]):

a=exp(tcoeff(%)), b=coeff(%,y1), c=coeff(%,z1);

            a = 0.0004189912543,  b =0 .996002908082386118, c = 0.0849219222885009006

L1:=x + y - 1 = 0, x - z - 1 = 0:

L2:=x + y - 7 = 0, x - y + 1 = 0:

A:=plots[intersectplot](L1, x=-5..5, y=-5..5, z=-3..5, thickness=2, color=red):

B:=plots[intersectplot](L2, x=-5..5, y=-5..5, z=-3..5, thickness=2, color=blue):

plots[display](A, B, scaling=constrained, axes=normal);