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.

The command zip is like map in that it provides elementwise operation. It differs from map in that it maps a two-argument function over a pair of structures. LinearAlgebra:-Zip is just like zip, but it also has an inplace option. So, you can achieve what you want with the single command

LinearAlgebra:-Zip(eval, BB, map[2,inplace](`=`, theta, CC), inplace);

Note that I used eval instead of subs, so the above produces evaluated results. If you actually want unevaluated results (perhaps to save processor time), do

LinearAlgebra:-Zip(subs, map[2,inplace](`=`, theta, CC), BB, inplace);

Why do you need to define f at all? You can use any symbol as a function without defining it. Just type f(y). But if you have a good reason for defining it, you can do

f:= x-> 'procname'(args);

If you want to insist that f only take exactly one argument, then do

f:= proc(x,$) 'procname'(x) end proc;

Now an error message will be issued if f is invoked with other than one argument.

plots:-spacecurve([cos(t),sin(t),t], t= 0..6*Pi);

The for loop

for i from a by b to c while d do
     e
end do

is equivalent to the do loop

i:= a:
do
     if (i > c) or not d then  break  end if;
     e;
     i:= i+b
end do;

The problem is that your function is undefined at 0. Are you looking for a rational function that closely approximates M6? We can expand at a=1 rather than a=0:

numapprox:-pade(M6, a=1, [3,3]);

The [3,3] indicates the degree of the numerator and denominator respectively.

There are several cases to consider:

The number of indices is known at the time the program is written: In this case, use nested loops:

for i from 1 to Ni do
     for j from 1 to Nj do
          ...
             for n from 1 to Nn do
                  ...
             end do
        ...
     end do
end do;

More generally, if the index values are stored in sets or lists  Si, Sj, ..., Sn:

for i in Si do
     for j in Sj do
         ...
               for n in Sn do
                       ...
                   end do
            ... 
    end do
end do;

The next two cases involve Cartesian products of the sets or lists of indices.

The number of indices is unknown at the time the program is written: In this case, use an iterator over the Cartesian product:

Iter:= combinat:-cartprod([Si, Sj, ..., Sn]);
while not Iter[finished] do
     v:= Iter[nextvalue]();
     # Now use v[1] as i, v[2] as j, etc.
     ...
end do;

The number of indices is unknown at the time the program is written and the Cartesian product is relatively small (say, less than 10 million entries): In this case, we generate the entire Cartesian product at once, because it is faster to do so:

CartProdSeq := proc(L::seq({set,list}))
local Seq,i,j;
option `Copyright (C) 2007, Joseph Riel. All rights reserved.`;
    eval([subs(Seq=seq, foldl(Seq
                              , [cat(i,1..nargs)]
                              , seq(cat(i,j)=L[j],j=nargs..1,-1)
                             ))]);
end proc:

for v in CartProdSeq(Si, Sj, ..., Sn) do

    # Now use v[1] as i, v[2] as j, etc.
     ...
end do;

The flush: The fancifully named pigeonhole principle states that if you put p pigeons into h holes then at least one hole will have at least x = ceil(p/h) pigeons. The 17 cards are the pigeons and the 4 suits are the holes: ceil(17/4) = 5. So how can we compute p = 17 knowing x = 5 and h = 4? p = (x-1)*h + 1.

The probability of a flush in a five-card draw is not as you state. It is binomial(4,1)*binomial(13,5)/binomiql(52,5). Permutations are not used because order does not matter.

The straight: Suppose that you had all the cards excepts the 5s and 10s. Then you'd have 44 cards and no straight. Add one 5 or 10 and now you have 45 cards and a straight (5*4^4 = 1280 straights in fact).

The probability of a straight in a five-card draw is 10*4^5/binomial(52,5).

You need to replace print(y) with print(eval(y)).

I don't think that I can give an adequate explanation of why. But try reading ?eval, beginning with the seventh paragraph under Description, the one beginning "The two remaining calling sequences...." The paragraphs above that apply to the two-argument version of eval, which is an unrelated command.