plot3d({8 - x^2 - y^2,x^2 + y^2}, x=-2..2, y=-sqrt(4-x^2)..sqrt(4-x^2));

The plot shoul look like

obtained with

plots:-inequal([x^2+y^2 >= x, x^2+y^2 <= 1/9], x=-1/3..1/3, y=-1/3..1/3,
optionsfeasible=[color=red], optionsclosed=[color=white], scaling = constrained);

The option coords=...  fails for some plot commands.

It's a version of the classical:

MDS:=proc(A::Matrix, r::integer) 
local n,B,i,j,ip,jp,k,P,W, rng;
rng:=t -> (t-1)*r+1 .. t*r; 
n:=upperbound(A,1)/r;
for i to n do   for j to n do 
  B[i,j]:=LinearAlgebra:-SubMatrix(A, rng(i), rng(j), datatype=integer[8])
od od;
for k to n do
   W := Matrix(r*k, datatype=integer[8]); 
   P:=combinat:-choose(n,k);
   for ip in P do for jp in P do
      for i to k do for j to k do
         W[rng(i),rng(j)] := B[ip[i],jp[j]];
      od od;
      if LinearAlgebra:-Modular:-Determinant(2,W)=0 then return 'NoMDS' fi;
   od od;
od;
'IsMDS'
end proc:

A:=Matrix(   # check
 < 1,0,1,0;
   0,1,0,1;
   1,0,0,1;
   1,1,1,0>):
MDS(A,2);
A[1,1]:=0:
MDS(A,2);
                             IsMDS
                             NoMDS

If you are working in Z2 (i.e. mod 2) then a square matrix A is MDS iff A = [1]  (so the dimension is 1).
MDS matrices are useful only for fields other than GF(2).
So, the MDS procedure reduces to

MDS:=proc(A::Matrix, r::integer)
  if r<> 1 then return 'NONE' fi;
  eval(A, 0=NO)
end proc;

e_n_1b := n_1 = (-w^(2*sigma)*tau + w^sigma)*s*nsp_1/(w^sigma*tau^2 - w^(2*sigma)*tau + w^sigma - tau) + tau*(w^(sigma - 1)*tau - w^(2*sigma - 1))*s*nsp_2/(w^sigma*tau^2 - w^(2*sigma)*tau + w^sigma - tau):
lhs(e_n_1b)=map(factor@expand, rhs(e_n_1b))

You want to simplify a polynomial f in n (=7) indeterminates.
Maple has simplify(f, size)  for this (size is by default).
Of course the size is not necessarily minimal, and AFAIK such algorithm does not exist.
But if you know/hope/suspect a special convenient form for f, tell us about it.

restart;
f:=
arccos((-sin(theta3)*cos(phi3)*cos(theta2)+cos(theta3)*sin(theta2))/((
(1-cos(theta3)^2)*cos(theta2)^2+cos(theta3)^2-1)*cos(phi3)^2-2*sin(theta3)*cos(
theta2)*cos(theta3)*sin(theta2)*cos(phi3)-cos(theta2)^2*cos(theta3)^2+1)^(1/2))
*sin(theta3):

evalf[15](Int(Re(f), [theta2=0..Pi,theta3=0..Pi,phi3=0..2*Pi], method=_CubaCuhre, epsilon=1e-6));
                        62.0124733278130
evalf[15](Int(Im(f), [theta2=0..Pi,theta3=0..Pi,phi3=0..2*Pi], method=_CubaCuhre, epsilon=1e-6));
                       Float(undefined)

See here.
It will be easy to use Explore and set your prefered sliders.

It is a bug in the symbolic part. The integral is simplified to BesselI for noninteger indices.

restart;
F := Int(exp(z*cos(theta))*cos(v*theta)/Pi, theta = 0 .. Pi):
ex := eval(F, [z=2, v=6/5]):
value(ex);     # wrong
    BesselI(6/5, 2)
evalf(%) = evalf(ex);
                   1.389981292 = 1.377481960
# FunctionAdvisor(BesselI, integral_form);

P.S. If you are intersted only in numeric results, just use Int instead of int and it will be OK.

It is used the fact that A,B are similar iff the characteristic matices (A-x*I, B-x*I) have the same Smith form.
So, the eigenvalues are not needed, and the computations are in Q if A,B are over Q.

A:=n ->
Matrix( n+1, (i,j)->
`if`(i<j,0,
  `if`(i=1,1,
  (i-1)*(-1)^(i-j)*2^(2*j-2)*(i+j-3)! / ( (2*j-2)!*(i-j)! )
))):

A(5);

You must use

PA[i] := eval(x[i + 1](t), sol(t[i]));

(note the i+1 because  x[0] is not defined) and similar for the rest of the loops.

Also, dt := 0.002;  because for 0.02 you will have singularities.

It is recommended to change tt[i] instead of t[i], because t is already in use.