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

restart;
evalc(exp(I*x));

                                           cos(x)+I*sin(x)

Or a prettier view:

restart;
with(InertForm):
a := 3;
b := 5;
c := -1;
eq := Display(exp(`%-`(c))/(a%*b) + a%^b*exp(c),inert=false);

This is similar to the classification of equilibrium points of two-dimensional linear dynamical systems. Below is an example of plotting several integral curves for the case of a stable focus:

restart;
DESys:={diff(x(t),t)=-x(t)-4*y(t), diff(y(t),t)=4*x(t)-2*y(t)};
DETools:-DEplot(Sys,{x(t),y(t)}, t=0..5, x = -1 .. 1, y = -1 .. 1, [[x(0) = 1, y(0) = 0],[x(0) = 0, y(0) = 1],[x(0) = -1, y(0) = 0],[x(0) = 0, y(0) = -1]]); 

There is no need to write the procedure choosing. Your example can be implemented in Maple as

L:=`$`~([w,b,g],4);
combinat:-choose(L, 4);

                L := [w, w, w, w, b, b, b, b, g, g, g, g]
 [[b, b, b, b], [b, b, b, g], [b, b, b, w], [b, b, g, g], [b, b, g, w], [b, b, w, w], [b, g, g, g], [b, g, g, w], [b, g, w, w], [b, w, w, w], [g, g, g, g], [g, g, g, w], [g, g, w, w], [g, w, w, w], [w, w, w, w]]

To solve the second problem, use the  ListTools:-Occurrences  command, for example:

e:=[b, b, b, w]: L0:=[w,b,g]:
[seq(s=ListTools:-Occurrences(s,e), s=L0)];

                             [w = 1, b = 3, g = 0]

Edit.

The procedure  Ellipse  does what you want:

restart;
Ellipse:=proc(a,b,phi)
local Vdot, dist, ell, F1, F2, M, MF2, MF1, varphi, c, m, t0, ELL, tp, po;
uses plots;
c := sqrt(a^2-b^2);
m:=2*arctan(c/b);
print(`maxϕ`=m);
if is(phi>m) then error "Should be phi<=`max&varphi;`" fi;
Vdot := proc(U, V) add(U[i]*V[i], i = 1 .. 2) end proc:
dist := proc(M, N) sqrt(Vdot(M-N, M-N)) end proc:
ell := x^2/a^2+y^2/b^2 = 1: 
c := sqrt(a^2-b^2):
F1 := [c, 0]: F2 := [-c, 0]:
M := [a*cos(t), b*sin(t)]:
MF2 := dist(M, F2);
MF1 := dist(M, F1);
varphi := arccos((MF1^2+MF2^2-4*c^2)/(2*MF1*MF2));
t0:=fsolve(varphi=phi, t=0..Pi/2);
print(parse("M")=eval(M,t=t0));
print(parse("MF2")=evalf(eval(MF2,t=t0)));
print(parse("MF1")=evalf(eval(MF1,t=t0)));
ELL := implicitplot(ell, x = -a-2 .. a+2, y = -2-b .. b+2, color = blue):
tp := textplot([[F1[], "F1"],[F2[],"F2"],[eval(M,t=t0)[],"M"]], 'align' = 'above', font=[Times,15]):
po := plot([F1,F2, eval(M,t=t0)], style = point, symbolsize = 12, symbol = solidcircle, color = red):
display([ELL, tp, po], scaling = constrained, size=[800,400]);
end proc:

Example of use:

Ellipse(7,5,Pi/4);

Edit. The procedure  Ellipse  text has been updated. I have added the calculation of the maximum angle  `max&varphi;`  and an error message if  phi > `max&varphi;`