If you don't mind extra parentheses in the display, you can do this

``(3) + ``(7) + ``(1);

If you want it without the parentheses, then it is significantly more complicated, but it can be done with the Typesetting package.

Your equation is fifth order in alpha. There is provably no possible general solution for fifth-order and higher-order polynomials. Since you have parameters, there is no numerical solution either.

When Maple makes up a variable name, it begins that name with a underscore to help distinguish it from user-supplied names. _Z takes the place of the solved-for variable in a RootOf expression. That is, _Z takes the place of alpha. A RootOf is Maple's placeholder for the root of an equation that Maple in unable to solve.

First, you need to change your syntax. To apply sqrt to two arguments, the closest syntax to what you had is

sqrt~([c^2, r^2]);

The reason that Maple does not simplify these square roots is that it does not know the sign of c and r or even that they are real. You can supply this information with an assuming clause:

sqrt~([c^2, r^2]) assuming c > 0, r > 0;

or

sqrt~([c^2, r^2]) assuming positive;

There is also a symbolic argument that can be give to sqrt:

sqrt~([c^2,r^2], symbolic);

add(a[i]*mul(`if`(j=i, 1, s-b[j]), j= 1..n), i= 1..n);

To obtain better understanding of what minimize does, enter your expression in exact form:

f:= x*sin(4*x)+11/10*sin(2*y):minimize(f, x= 0..10, y= -10..10, location);

Note that the returned answers are exact. I think that it finds the zeros of the gradient and evaluates the function at them to find the minimum evaluation. It does not use a numerical method.

Answering your second question, Explore can be used with piecewise just as with any other mathematical expression. For example,

Explore(plot(piecewise(x < a, -1, 1), x= -10..10, discont), parameters= [a = -10..10]);

I don't know what's going with your first question. It always appears in the same document for me. Maybe there is a default setting for it.

The command with does not work correctly inside procedures. Eliminate it, and replace the jacobi line with

numtheory:-jacobi(A,B)

Indeed, the entire procedure is superfluous, and could be replaced by the line

jacob:= numtheory:-jacobi;

First, remove Pr from your list of parameters in the second line. Second, change the line where you define PDE to

PDE:= unapply({pde||(1..3)}, Pr);

Then replace your last two lines with

PrList:= [0.71, 7, 10]:
Colours:= table(PrList =~ ["blue", "red", "green"]):
for Pr in PrList do
     pds:= pdsolve(PDE(Pr), IBC, numeric);
     Plots[Pr]:= pds:-plot[display](
          u(y,t), t=0.5, linestyle= "solid", colour= Colours[Pr],
          legend= sprintf("Pr=%2.2f", Pr),
          title= "Velocity Profile", labels= ["y", "theta"]
     )
end do:
plots:-display([seq(Plots[Pr], Pr in PrList)]);

There is something weird going on when Gamma is an odd integer, including Gamma=1.

To explore further, first let's remove the extraneous and obsfucating elements. This was afterall a coloring function that I wrote for a highly specialized situation. Here's a simplified version animated over a range of values of Gamma.

restart:
plots:-display(
     [seq(
          plot3d(
               0, 0..1, 0..1,
               color= [
                    (x,y)-> (1-x)^Gamma/3, #Hue
                    (x,y)-> 1,             #Saturation
                    (x,y)-> 1,             #Value
                    colortype= HSV
               ],
               title= sprintf("Gamma = %1.1f", Gamma)
           ),
           Gamma= 1..6, .2
      )],
      insequence,
      lightmodel= NONE,
      style= patchnogrid,
      orientation= [-90,0]
);

I recommend single-frame-stepping through the animation. Notice what happens at odd integer values of Gamma. The denominator 3 in hue is supposed to select for the lower third of the visual spectrum. And it works fine except for those odd integer values of Gamma for which the whole spectrum is displayed. Changing the 3 to any other value greater than 1 doesn't help. Forcing Gamma to floating point doesn't help. Changing 1-x to x does correct the situation (and, of course, reverses the direction of the colors).

I don't know why. Acer, do you have any idea? You may be the only person here who can figure this out.

I changed capital I to capital J because I'd rather just avoid using capital I as a variable.  I made the main loop start at 0 instead of 1. The 0 didn't make sense to me.

restart:
T[0]:= 5.5556e7:  J[0]:= 1.1111e7:  V[0]:= 6.3096e9:
A1:= 1:  A2:= 1:  
c:= 0.67:  h:= 1:  d:= 3.7877e-3:  delta:= 3.259*d:
lambda:= 2e8/3*d:  R0:= 1.33:
p:= c*V[0]*delta*R0/lambda*(R0-1):
beta:= d*delta*c*R0/lambda*p:
#Step 1
n:= 100:
lambda__1[n]:= 0:  lambda__2[n]:= 0:  lambda__3[n]:= 0:
u__1[0]:= 0:  u__2[0]:= 0:
#Step 2:
for i from 0 to n-1 do
     T[i+1]:= (T[i]+h*lambda)/(1+h*(d+(1-u__1[i])*beta*V[i]));
     J[i+1]:= (J[i]+h*(1-u__1[i])*beta*V[i]*T[i+1])/(1+h*delta);
     V[i+1]:= (V[i]+h*(1-u__2[i])*p*J[i+1])/(1+h*c);
     lambda__1[n-i-1]:= (lambda__1[n-i]+h*(1+(1-u__1[i])*beta*V[i+1]))/
          (1+h*(d+(1-u__1[i])*beta*V[i+1]));
     lambda__2[n-i-1]:= (lambda__2[n-i]+h*lambda__3[n-i]*(1-u__2[i])*p)/(1+h*delta);
     lambda__3[n-i-1]:= (lambda__3[n-i]+h*(lambda__2[n-i-1]-lambda__1[n-i-1])
          *(1-u__1[i])*beta*T[i+1])/(1+h*c);
     R1:= (1/A1)*(lambda__1[n-i-1]-lambda__2[n-i-1])*beta*V[i+1]*T[i+1];
     R2:= -(1/A2)*lambda__3[n-i-1]*p*J[i+1];
     u__1[i+1]:= min(1, max(R1,0));
     u__2[i+1]:= min(1, max(R2,0));
end do;

restart:
X:= Statistics:-Distribution(ProbabilityFunction= (x-> (2*x+1)/25), Support= 0..4):
Statistics:-CDF(X, x);

A more primitive method:

f:= x-> (2*x+1)/25:
sum(f(k), k= 0..x);

1. You say "I need to call a[j], b[j], E[j][1], E[j][2], and the coordinate points." What are the coordinate points? Is there one for each row?

2. The last three points in array P are never used. Is that correct?

3. Each loop iteration defines only one of a[j] and b[j]. What to do about the other?

Here's how you can store the computed values in a Matrix. I just added two lines to your code, which I've highlighted below:

P:= [[8, 4], [8, 3], [8, 2], [7, 1], [6, 0], [5, 0], [4, 0], [2, 1], [1, 1], [1, 4]];
M:= Matrix(4,5);
for j from 2 to 5 do
     k[j] := j+1;
     x[j] := add(P[j, 1], j = j-1 .. j+2);
     X[j] := add(P[j, 1]^2, j = j-1 .. j+2);
     y[j] := add(P[j, 2], j = j-1 .. j+2);
     Y[j] := add(P[j, 2]^2, j = j-1 .. j+2);
     xy[j] := add(P[j, 1]*P[j, 2], j = j-1 .. j+2);
     cx[j] := evalf(x[j]/k[j]);
     cy[j] := evalf(y[j]/k[j]);
     c11[j] := evalf(X[j]/k[j]-cx[j]^2);
     c22[j] := evalf(Y[j]/k[j]-cy[j]^2);
     c12[j] := evalf(xy[j]/k[j]-cx[j]*cy[j]);
     C[j] := evalf(Matrix(2, 2, [[c11[j], c12[j]], [c12[j], c22[j]]]));
     E[j] := simplify(fnormal(LinearAlgebra[Eigenvalues](C[j])));
     if E[j][1] > E[j][2] then
          a[j]:= E[j][2]/(E[j][1]+E[j][2])
     else
          b[j]:= E[j][1]/(E[j][1]+E[j][2])
     end if;
     M[j-1, 1..4]:= < a[j], b[j], E[j][1], E[j][2] >
 end do;

Even though it doesn't cause any problems in your code, it is generally considered a bad practice to use the same variable as both the index in an add and in the bounds of the add. But I did not change your code in the several places that you did this.

Here are the very basics of logarithmic rescaling to get an even grid spacing: Let's say that your original plot is

densityplot(f(x,y), x= a..b, y= c..d);

and you want logarithmic scaling on the first axis. Then change the plot command to

densityplot(f(10^x, y), x= log10(a)..log10(b), y= c..d);

If f is a procedure, then it is trickier, but still based on the same idea.

Then you need to redo the tickmarks.

I made some corrections to your worksheet, as noted by the comments in it. I am sure about most of the corrections. I also removed the principal normal computation---I am not sure about that. I think that PathInt takes it into account automatically. If you know what final results you expect, then you can verify with what I got. If you still need the principal normal , I can work on that.

Your worksheet was in fact using exact arithmetic, and I corrected that.

(Partially?) corrected worksheet: main_screened_Poiss.mw

The problem is the quote marks that you use on ':-language'= language. I have used the correct quote marks here. The quotes that you used do not occur on the American keyboard, and I can't reproduce them here. Needless to say, they are not part of the Maple language. The correct quote mark is character number 39 in ASCII.