I removed unnecessary options:

restart;
F1:=1/(x + mu) + x - 1;
plots:-implicitplot(F1 = 0, mu = -8 .. 8, x = -4 .. 4, color=blue, thickness=2, gridrefine = 3);

Let a point move in the positive direction along the unit circle. It makes a complete revolution in time  T seconds and stops. Let its angular velocity decreases uniformly from the value  v0  to 0. Below we find the law of motion (the dependence of the angle of rotation  omega on the time  t) by solving the differential equation  Eq  and show the  N  positions of this point at regular time intervals and its radius-vector. Animation of this movement is also made. Of course we can change the parameters  T  and  N  as we like.

restart;
T:=10:
v:=t->-k*t+v0: # Angular velocity as a function of time t
# Law of motion: omega(t) - angle of rotation as a function of time
Eq:=diff(omega(t),t)=v(t);
dsolve({Eq,omega(0)=0});
omega:=unapply(eval(omega(t),%), t);
Sys:={v(T)=0,omega(T)=2*Pi};
solve(Sys,{v0,k});
assign(%):

with(plots): with(plottools):
N:=20:
Circle:=circle(color=blue):
t:=T*i/N:
Radius_Vector:=seq(line([0,0],[cos(omega(t)),sin(omega(t))],color=red,thickness=2), i=0..N-1):
Point:=seq(disk([cos(omega(t)),sin(omega(t))],0.025,color=red), i=0..N-1):
display(Point,Circle,Radius_Vector, axes=none);

P:=proc(t)
uses plots, plottools;
display(line([0,0],[cos(omega(t)),sin(omega(t))],color=red,thickness=2),disk([cos(omega(t)),sin(omega(t))],0.025,color=red));
end proc:

animate(P,['t'], 't'=0..10, frames=100, background=display(Circle), scaling=constrained, trace=20
);

radius_vector.mw

See help on  timelimit  command.

Denoting the elements of the matrix  T__PhPh  by  x[i,j] , i=1..3, j=1..3, we obtain a system of 3 linear equations with 9 unknowns. Obviously, in the general case, such a system has an infinite number of solutions depending on 6 parameters.

Q20210127_new.mw

restart;
M:=Matrix(3, 3, [[-1/3*a + b + c - 2*d, -1/4*a + 2/3*c, 1/15*a - 1/4*c + 2/3*d], [e + a - 2*c, 2/3*a, -1/4*a + 2/3*c], [f - 2*a, 0, 2/3*a]]);
Id:=Matrix(3,3,(i,j)->`if`(i=j,1,0));
Sys:=[seq(seq(M[i,j]=Id[i,j],j=1..3),i=1..3)];
solve(Sys);

Similarly, you can solve for any  n .

remove(t->is(op(1, t) = -1), [T]);

See help on the command  Student:-Calculus1:-ApproximateInt  for this.

Before looking for a root, it is always useful to do the simplest qualitative analysis of the equation. Since we are looking for real roots, we first find the domain of the expression  P . Obviously should be  6*lambda - 2 >= 0  or  lambda>=1/3 . So it doesn't make sense to consider  lambda < 1/3  where the expression  P  takes complex values (this is the reason for the error in your file) .

Having built the plot, we clearly see one root. This root is easily found by  fsolve  command that Carl did. Looking at the plot it seems that this root is the only one, but of course this requires a separate proof, since any plot shows the behavior of a function only in a certain finite interval. In this example, one can prove the uniqueness of the root simply by finding the derivative of  P  and checking that it is negative for all  lambda>1 . Below we also show 3 more ways to find this root.

restart;
P := (-216*sqrt(2)*lambda^2*(lambda-4/3)*arctan((1/2)*sqrt(6*lambda-2)*sqrt(2))-108*Pi*lambda^2*(lambda-4/3)*sqrt(2)+(-12*lambda+16)*(6*lambda-2)^(3/2)+(-40*lambda+32)*sqrt(6*lambda-2))/(576*lambda^2)+1/(6*lambda);
plot(P, lambda=0..2, size=[600,250]);
simplify(diff(P, lambda)); # The derivative of P
r:=fsolve(P=0); # The unique real root

# Other ways to find the root
RootFinding:-Analytic(P, re = 0.5 .. 10, im = -1 .. 1); 
Student:-Calculus1:-Roots(P, lambda = 0.5 .. 10, numeric)[]; 
Student:-NumericalAnalysis:-Bisection(P, lambda = [0.5, 10], tolerance = 10^(-9));

restart;
A:=[a,-b,c, d, -a,b, -d];
remove(t->sign(t)=-1, A);
convert(%, set);

                                 A := [a, -b, c, d, -a, b, -d]
                                            [a, c, d, b]
                                            {a, b, c, d}

restart;
A := [x[1, 1]*(a*x[1, 1] + E[1, 1]) + x[2, 1]*(a*x[1, 2] + E[1, 2]) + x[3, 1]*(a*x[1, 3] + E[1, 3]), x[1, 2]*(a*x[1, 1] + E[1, 1]) + x[2, 2]*(a*x[1, 2] + E[1, 2]) + x[3, 2]*(a*x[1, 3] + E[1, 3]), x[1, 3]*(a*x[1, 1] + E[1, 1]) + x[2, 3]*(a*x[1, 2] + E[1, 2]) + x[3, 3]*(a*x[1, 3] + E[1, 3]), x[1, 1]*(a*x[2, 1] + E[2, 1]) + x[2, 1]*(a*x[2, 2] + E[2, 2]) + x[3, 1]*(a*x[2, 3] + E[2, 3]), x[1, 2]*(a*x[2, 1] + E[2, 1]) + x[2, 2]*(a*x[2, 2] + E[2, 2]) + x[3, 2]*(a*x[2, 3] + E[2, 3]), x[1, 3]*(a*x[2, 1] + E[2, 1]) + x[2, 3]*(a*x[2, 2] + E[2, 2]) + x[3, 3]*(a*x[2, 3] + E[2, 3]), x[1, 1]*(a*x[3, 1] + E[3, 1]) + x[2, 1]*(a*x[3, 2] + E[3, 2]) + x[3, 1]*(a*x[3, 3] + E[3, 3]), x[1, 2]*(a*x[3, 1] + E[3, 1]) + x[2, 2]*(a*x[3, 2] + E[3, 2]) + x[3, 2]*(a*x[3, 3] + E[3, 3]), x[1, 3]*(a*x[3, 1] + E[3, 1]) + x[2, 3]*(a*x[3, 2] + E[3, 2]) + x[3, 3]*(a*x[3, 3] + E[3, 3])];

subs({seq(E[k,k]=1, k=1..3)}, A);

Below is another way that more closely matches your request:

subsindets(A, Or(seq(identical(a*x[k,k]+E[k,k]), k=1..3)), t->op(1,t)+1); 

Your equation  pg1  contains 3 parameters  k, Q, lambda  besides the unknown sigma. Indicate the values of these parameters and solve numerically by  fsolve  command. For example:

fsolve(eval(pg1, [k = 2, Q = 3, lambda = 4]) = 0, sigma);

                                       -11.63702103

restart; 
rho := 1/2: mu := 1/2: T := (1/4)*(t-x)^(1/2): 
E := add(T^k/GAMMA(k*rho+mu), k = 0 .. 5); 
eq1 := h(x) = 8/3*(int((t-x)^(1/2)*E*(1/(sqrt(Pi)*x^(1/2))-3/32*(x^2)+3/16*exp(-x)-1), t = -1 .. 1)): 
eq2 := simplify(eq1);

Edit. Your expression  h(x)  takes complex values for any real values of  x . Therefore, you cannot plot  y=h(x)  for real  x  from a certain interval, but you can plot, for example,  y=abs(h(x)) :

plot(abs(rhs(eq2)), x=-3..3, 0..10);

It happens that Maple, for some calculations, returns the result as a float-number with a very small imaginary part. If you are sure that really the result is a real number, then you can get rid of the imaginary part by using the commands  fnormal  and  simplify  with the  zero  option (see help on these commands for details).

Example:

1.5+2.e-15*I;
fnormal(%);
simplify(%,zero);

The  D  is the differentiation operator (see help for detailes). In my opinion, the following form gives a better understanding of the  Chain rule :

restart;
diff(V(h(x)),x);
convert(%, diff);

Use more intelligent  &^  instead of  ^  for this:

restart;
n:=4142074788597595058238933087635402958555715420375128103710780511:
m:=65108105099101044032:
e:=3467995563410655813589059428902496618444213189544232426283402949:
m &^e mod n;

                 24733671999927750262610339317849962811213365592569778941953149