Obviously the OP Rodey Genius is not interested in learning. However, I had this Answer already 80% done at the time that they deleted the Question, so I'd rather not waste what I coded. Maybe someone else could learn from it.

I chose to code determinants by cofactor expansion because it can be used to illustrate a great many features of Maple programming in one short procedure. The algorithm is straightforward to code recursively. I optimized it the way one would do by hand for a sparse matrix: At each iteration, I expand along the row or column with the most 0s.

I think that even most Maple experts will learn something by studying this procedure.

As a bonus, I found an (admittedly, rather obscure) class of matrices for which my procedure significantly outperforms default methods: half-sparse matrices of bivariate rational functions.

Determinant:= proc(
    A::Matrix(square), 
    {method::identical(cofactors, other):= 'other'}
)
uses AT= ArrayTools, St= Statistics, LA= LinearAlgebra;
local k, n, r, c, T, transp, rm1, rp1, a;
    if method='other' then return LA:-Determinant(A) fi;
    
    # method = cofactors:
    n:= upperbound(A)[1];
    if n<2 then return `if`(n=0, 1, A[1,1]) fi; #base cases
    #Count number of non-zero entries in every row and column:
    T:= St:-Tally~([((r,c):= AT:-SearchArray(A))]);
    if nops~(T) <> [n,n] then return 0 fi; #There's a row or column of 0s.
    T:= sort~(T, 'key'= rhs);
    transp:= evalb(rhs(T[1][1]) > rhs(T[2][1])); #expand a row, or a column?
    r:= lhs(T[`if`(transp, 2, 1)][1]);  rm1:= r-1;  rp1:= r+1;
    add(
        if (a:= `if`(transp, A[k,r], A[r,k])) = 0 then 0
        else
            (-1)^(r+k)*a*thisproc(
                A[[[..rm1,rp1..],[..k-1,k+1..]][`if`(transp, [2,1], [1,2])][]], 
                _options
            )
        fi, 
        k= 1..n
    )
end proc
:
#Test case: half-sparse matrix of bivariate rational functions.
A:= LinearAlgebra:-RandomMatrix(
   9, density= 0.5, 
   generator= (()-> randpoly(x, degree=1)/randpoly([x,y], degree= 1))
):
d1:= CodeTools:-Usage(Determinant(A, method= other)):
memory used=23.74MiB, alloc change=0 bytes, cpu time=938.00ms,
real time=482.00ms, gc time=0ns
d2:= CodeTools:-Usage(Determinant(A, method= cofactors)):
memory used=11.99MiB, alloc change=0 bytes, cpu time=109.00ms, 
real time=130.00ms, gc time=0ns
expand(numer(d1-d2)); #Verify equality of results.
                               0
#Verify results are not identically 0:
eval(d2, [x= 0, y= 0]);
           16003080668688229738162772913864379992367
           -----------------------------------------
            86804823806332323196340023957536768000  

There's nothing wrong with Kitonum's Answer, but this alternative form is often useful, especially when n is not integer: For x, y, and n real,

(x+I*y)^n = (x^2+y^2)^(n/2)*exp(I*n*arctan(y,x))

Easy analytic geometry solution:

A:= <1,2,3>:  v:= <1,2,2>:  C:= <10,11,12>:
L:= A+t*v:

We want the point where line L intersects the plane normal to L containing point C:

B:= eval(L, t= solve(L.v = C.v));

Plot:

plots:-display(
    plots:-implicitplot3d(
        <x,y,z>.v = C.v, ([x,y,z]=~ 0..15)[], 
        style= surface, transparency= .2
    ), 
    plots:-pointplot3d([A,B,C], connect, color= red, thickness= 3),
    orientation= [80, 80, 130], axes= frame, labels= [``$3]
);

Download BasicLA.mw

If you use option matrixview, then the vertical-axis tickmarks are indexed -1, -2, ..., -20.

The command SPolynomial is in package Groebner not package Ore_algebra.

You need to do two things:

1. Your plot range 0..5 is much too big for this problem. Change it to approximately 0..0.01 (for all plots).

2. Add the option maxfun= 0 to the dsolve command. This will remove any limitation on the number of points that it needs to compute.

3. (Optional) Remove output= listprocedure. It doesn't do anything beneficial in this case, and I'm always suspicious of it.

Here it is. Obviously the placement of the sliders could be prettier, but I think mathematically it's all done. MaplePrimes can't handle Explore plots, but it should be there when you download the worksheet.
 

Download VaccinationRateSIR.mw

You're missing a closing parenthesis ) at the end of the line immediately above end. Note that when your cursor is immediately after a closing bracket, then a faint box appears around it's opening match. 

Looks like you're doing the elliptic curve addition! Good.

You have a few more errors that will just cause weirdness rather than an error message:

1. Assignment is := not =.
2. Expand does nothing without mod.
3. Multiplication should use an *.

Here is a phase portrait for the three-tree-frog system using some parameter values and initial conditions that I just made up. You'll need to do a lot of playing around with the parameters to figure out the rest of part (c) of the problem.
 

Download ThreeTreeFrogs.mw

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);