Tima, I don't think that there is a simpler way to solve this, but there are a few much better ways. The problem with most standard solutions, such as the one you give, is that the real roots end up with false small imaginary parts and it's iffy whether those roots get plotted. (I call them spurious imaginary parts.) Here are two methods that avoid that problem. The first method also gives exact symbolic solutions (just as yours does), but they are constructed such that evaluating them numerically gives pure real numbers.

To understand the rest of this, it helps to know that there are 2 real roots for a=0 and 3 real roots for any a in (0,1]. It's easy to prove this with elementary calculus, and it's also verified by solve(..., real, parametric) solution. You should look at that solution in the worksheet, but it's so wide on the screen that I don't want to put it directly on MaplePrimes.

Download SolveParametric.mw

restart:
f:= (x,a)-> x^3 + (a-3)^3*x^2 - a^2*x + a^3:
Titles:= ["Least root", "Middle root", "Greatest root"]:
PlotIt:= R-> 
    local k,
    Views:= `[]`~(`..`~(R~([0,1])[]), 0..1):
    plots:-display(
        `<|>`(
            seq(
                plot([a-> R(a)[k], a-> a, 0..1], view= Views[k], title= Titles[k]),
                k= 1..3
            )
        ),
        labels= [x,a], gridlines= false
    )
:
Method 1: solve(..., real, parametric): This gives the exact symbolic roots in 
RootOf form but in such a way that floating-point evaluation of the real roots 
does not produce spurious imaginary parts.
st:= time():
SolP:= solve(f(x,a) = 0, x, real, parametric):
SolPR:= subsindets(
    subs([x=0]= ([x=0]$2), SolP), #Account for double root at x=0, a=0.
    [identical(x)= anything],
    e-> op([1,2], e) #I.e., change every [x=r] to r.
):
P1:= PlotIt(proc(a) option remember; evalf(eval(SolPR, :-a= a)) end proc):
time() - st; 
                             2.922

Method 2: fsolve: This also avoids spurious imaginary parts (for 
polynomials with real coefficients). But it doesn't give exact, symbolic
solutions.
st:= time():
P2:= PlotIt(proc(a) option remember; local x; [fsolve(f(x,a))] end proc):
time() - st;
                             0.766

P1; #plot via solve(..., parametric)

P2; #plot via fsolve

J:= <
    118, 174, 374, 597, 670, 677, 516, 407, 
    291, 196, 155, 101, 66, 34, 21, 9
>:
A:= <seq(2.5+5*k, k= 0..15)>:
plot(
   <A|J>, style= point, symbol= asterisk, symbolsize= 16,
   labels= [Age, JOHOR], labeldirections= [horizontal,vertical],
   view= [0..80, 0..1.1*max(J)], axes= boxed
);

Mathematically speaking, the job can be done with a single short command, plottools:-transform. The shaded regions will turn out a bit ragged, but a mathematically sophisticated viewer will understand what's represented. An alternative is to not use shading in the original plot, and then the transform'd result will look perfect, like this:

restart:
z:= x+I*y: #Complex z must be defined in terms of real x and y for this to work.
A:= plots:-implicitplot(
   [evalc(abs(z)) < 3, evalc(1 < abs(z - 1))],
   x= -5..5, y= -5..5,
   scaling= constrained, grid= [100$2]
);
#Regardless of how you construct the original plot above, the code below
#remains the same:
x1:= (9 - sqrt(45))/2; x2:= (9 + sqrt(45))/2;
M:= (3 + sqrt(5))/2*(z - x1)/(z - x2);
plottools:-transform(unapply(evalc([Re,Im](M)), (x,y)))(A);

Note that the argument of transform is a function from R² to R² where the output is a list (i.e., in square brackets such as [x,y]) and the input is just a sequence (i.e., such as (x,y)). The output of transform is itself a function which can be applied to a plot---in this case, any 2D plot.

Using Maple 2019, if I give your second command, the one with discont, then I get exactly what you want. So, there's some bug in your Maple 2015 that has since been fixed. The same bug may also affect the output of the following plot command; on the other hand, it may work in Maple 2015:

plot(
   3/(1-3/2*sin(theta)), theta= 0..2*Pi, discont, 
   coords= polar, axiscoordinates= polar, 
   coordinateview= [0..10, 0..2*Pi], discont
);

This should produce a plot identical to the equivalent polarplot command, and indeed it does in Maple 2019.

All of the inputs to the GF field operations need to be processed by ConvertIn before use. (For what it's worth, this is needed to put them into the zppoly format that I mentioned earlier. But you don't need to understand that format to work with GF.) So, this works:

G:= GF(2, 3, T^3+T^2+1):
t:= G:-ConvertIn(T):
G:-`*`(t,t);

You also asked:

• how do I specify what variable is used for my finite field without needing to specify an irreducible polynomial?

The field elements are indexed in a predictable way. The index of the extension element itself is always the field characteristic (the prime). So, if you want to use T, do

restart:
G:= GF(2,3):
T:= G:-input(2):
G:-`*`(T,T);

One would think that if the "pattern" variables a, b, i were made local, that would solve the problem. It's easy enough to make them local:

EinsteinRule:= e->
    local a, b, i;
    applyrule(
        vec(a::symbol, i::symbol) * vec(b::symbol, i::symbol) = 
            dotproduct(a,b),
        e
    )
:

But I am absolutely dumb-founded that that doesn't fix your issue. There should now be absolutely no connection between the a in the applyrule command and any a that you may happen to use in the arguments. Yet...

EinsteinRule(vec(a,k)*vec(b,k));
                        dotproduct(a, b)
EinsteinRule(vec(u,k)*vec(u,i)*vec(c,k));
                   dotproduct(c, u) vec(u, i)
EinsteinRule(vec(a,k)*vec(a,i)*vec(c,k));
                 vec(a, k) vec(a, i) vec(c, k)

That's a critical flaw in applyrule. It's as if it's converting things to strings and comparing the strings.

I've never liked applyrule ​​​​​​anyway. Here's a version using the much more versatile command subsindets:

EinsteinRule:= e-> 
    subsindets(
        e,
        `*`,
        proc(p) 
        local V,v,o,pv,c;
            (V,o):= selectremove(type, p, :-vec(symbol$2));
            v:= {op(V)};
            for pv in combinat:-choose(v,2) do
                c:= op~(pv) minus `intersect`((`{}`@op)~(pv)[]);
                if nops(c)=2 then
                    V/= mul(pv)/dotproduct(c[]);
                    break
                fi
            od;
            V*o
        end proc
    )
:

Now, try the same three examples.

I am not familiar with the Einstein summation convention, so I tried to generalize reasonably just from the information that you provided. My procedure may not have captured every nuance of it.

If you read my Answer to the other similar Question that you Replied to, you'll see how to work with finite fields of non-prime order using just the mod command. This is by far the easiest way to start. After that, you can learn GF, then after that modp1 and modp2.

Put your cursor in the Code Edit Region and right-click to bring up the context menu. Select Execute. The results (if any, including error messages) will be displayed immediately after the Region. If there are no syntax errors, then you can now call your function from the worksheet.

The commands apply and unapply are red herrings here (idiom for "misleading clue"). They will not help you understand function application.

However, the help page ?evalapply will help a lot. But note that the (very old) command that it describes is almost entirely handled automatically in the kernel now. So, this page helps you understand what's happening rather than what you need to code to make it happen.

Here is the improved procedure. It's a little bit faster. Please test it. I don't know if it'll run in your version of Maple. I do think that it'll be necessary to run it in 1D-input (aka Maple Input).

BP:= proc(A::listlist({0,1}))
local
    n:= nops(A), m:= nops(A[1]), 
    S:= convert(Matrix((m,m), shape= identity), list, nested), 
    Q, R, T, j, QQ, Break, Z, Z2, E, uk:= 0, U:= Vector()
;
    for uk do
        E:= table();
        for Q in combinat:-choose(nops(S),2) do
            R:= [S[], `+`(S[Q][]) mod 2];
            (Z,Z2):= (0,0);
            for T in A do
                for j to nops(R) do 
                    for QQ in combinat:-choose(R,j) do                    
                        Break:= evalb(`+`(QQ[]) mod 2 = T);
                        if Break then  (Z,Z2):= (Z+j, Z2-j^2);  break  fi 
                    od;
                    if Break then break fi
                od
            od;
            E[Z, Z2, Q]:= ()     
        od;
        (Z,Z2,Q):= sort([indices(E)])[1][];
        S:= [S[], `+`(S[Q][]) mod 2];
        U(uk):= Q; 
        if Z=n and Z2=-n then return [seq(U)] fi
    od
 end proc
:

There are several ways that it could be done with indets. One way is

indets(P, function(Not(function)))

Specifically, this says to extract the functions none of whose arguments are themselves functions.

The problem with using select for this is that it only scans the top level of structures. A problem with using patmatch is that it tries to match the whole structure rather than subparts.

Anything at all can be made into a Maple symbol by enclosing it in back-quotes: `Diff(f(t),t)`.

The solve command can already solve for the highest-order derivatives without needing to use the above trick.

If the score is 0 for an unanswered question, then the possible scores are 

{seq(seq(3*c-i, i= 0..10-c), c= 0..10)};
{-10, -9, -8, -7, -6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 
8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23,  24, 26, 27, 30}
nops(%);
                               38

So 38 is the maximum number that could have distinct scores.

If, on the other hand, the score is -1 for both unanswered and incorrectly answered questions, then the possible scores are

{seq(3*c-(10-c), c= 0..10)};
          {-10, -6, -2, 2, 6, 10, 14, 18, 22, 26, 30}
nops(%);
                               11

So 11 is the maximum number that could have distinct scores.

If partial credits are given (1 or 2 points for partially correct answers), then of course there are 41 possible scores: -10, ..., 30.

fu_bar:= u-> u^2/2;
fu_barr:= fu_bar~(u0);

Using the arrow -> is equivalent to proc. The above is equivalent to

fu_bar:= proc(u) u^2/2 end proc:

No local variable is needed in this case, although what you had works anyway.

An option with no proc and no arrow, just the elementwise operator ~, is

fu_barr:= u0^~2/2;

It's usually better to perform operations on an entire matrix as a whole. In your case you could do

Am:= ArrayTools:-AddAlongDimension(A,1)/k:
C:= Matrix((k,k), (i,j)-> A[i,j] - Am[j]):