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ϕ`  and an error message if  phi > `maxϕ`

If I understood your question correctly, then you want to learn how to build such graphs in Maple. Your graphs are usual broken lines. The simplest way to plot them in Maple is simply to specify a list of vertices of the broken line. Below, as an example, in Maple, the first graph is built:

restart;
L:=[[0.5,0.43],[0.525,0.32],[0.55,0.28],[0.575,0.52]]:
plot(L, thickness=2, size=[500,300], title="Wnf", labels=["b",``], font=[times,16], labelfont=[times,18]);

restart;
Digits:=20:
L:=[x^2/5-y^2/2=1,y=x^3-7/2*x^2+2*x]:
f:=solve(L[1],x);
g:=solve(L[2],x);
Sol:=evalf(RealDomain:-solve(L));
plots:-implicitplot(L, x=-2..4,y=-3..3, color=[red,blue], gridrefine=3);
Area=simplify(fnormal(int(g[1]-f[1], y=op([2,2],Sol[2])..op([2,2],Sol[1]))), zero);

                                Area =0 .75833670466253534952

restart;
for i from 1 to 2 do
for j from 1 to 2 do
T[cat(i,j)]:=2*mu*varepsilon[cat(i,j)]-add(2*mu/3*varepsilon[cat(r,r)]*delta[cat(i,j)],r=1..3);
print(S[cat(i,j)]=T[cat(i,j)]);
od: od:

Or (all indices in roman):

restart;
f:=(x,y)->parse(convert(cat(x,y),string)):
for i from 1 to 2 do
for j from 1 to 2 do
T[f(i,j)]:=2*mu*varepsilon[f(i,j)]-add(2*mu/3*varepsilon[f(r,r)]*delta[f(i,j)],r=1..3);
print(S[f(i,j)]=T[f(i,j)]);
od: od:

restart;
ode:=x = (diff(y(x),x)^2+1)^(1/2)*diff(y(x),x);
mysol:=y(x)=_C1+int(  sqrt(-2+2*sqrt(4*a^2+1))/2,a=0..x);
check:=odetest(mysol,ode);
simplify(combine(check)) assuming x>0;