Your first plot specification is [X[k],Y[k]], k= 0..n. This type of specification treats k as a continuous variable. Since k is an integer variable, use seq, replacing the above with

[seq([X[k],Y[k]], k= 0..n)]

(I see now that the above is substantially the same as Kitonum's Answer.)

Here's another way, usings Arrays, and a single plot command:

f:= x-> x^2 + x - 12;
(a, n, h):= (0, 10, 5);
X:= Array(0..n, i-> a+(i-1)*h, datatype= hfloat); #Also works w/o datatype= hfloat.
Y:= map[evalhf](f,X); #Also works w/o [evalhf].
printf("\n    i        x       f (dec.form)    f (sci. notat.)\n"); 
printf( "  ---------------------------------------------------\n"); 
seq(printf(" %5d  %9.4f  %14.9f  %18.10e\n", i, X[i], Y[i], Y[i]), i= 0..n);
plot(
   # <...> is Matrix/Vector constructor/amalgamator.
   # (...)^+ is Matrix transpose.
   [<X,Y>^+, f(x)], x= a..a+n*h, 
   style= [point, line], symbol= soliddiamond, symbolsize= 24,
   color= [blue, black]
);

String literals, such as "A", should be in double quotes. Since you're also using A as a matrix name, it can't be a string also. But "A" is always a string.

You have two errors, each of which occurs twice. The first error is that if you're going to use both a "parent" variable (such as your f or theta) and a subscripted form of that variable (such as f[0] or theta[0]), then you can't make an assignment to either form or to a function made from it, such as f[0](eta). If you do it, it doesn't immediately cause an error; it just leads to erroneous results later^[*1]. So, in L[1] and L[2], you should change both f and theta to something else.

The second error is that in your assignments to ode1 and ode2, you need to put a space after L[1] and L[2] for there to be "implied multplication". Placing a name (such as L[1]) immediately before a left parenthesis is interpretted as function application rather than multiplication.

^[*1]In particular, making an assignment to f[0] turns f into a table, which is why you see table in the error message.

Your PDF file has one differential equation with two unknown functions: Z(t) and V(t). Thus, you'll need one more equation with one or both of those functions. You'll need two initial condition for Z (for example Z(0) and Z'(0)) and perhaps also one or more for V, depending on the highest differential order of V (if any) that occurs in the second equation.

If you have a sequence of equations EQs and a sequence of initial conditions ICs, then a plot of V vs. t can be obtained by simply:

sol:= dsolve({EQs, ICs}, numeric);
plots:-odeplot(sol, [t, V(t)], t= 0..2);

where you can change 0..2 to any range that you want.

The fact that the equations are nonlinear is irrelevant; the commands and the numeric solution techniques are the same regardless.

Suppose that the existing package (EP) is modified by its author (A). How will you handle that? Specifically, think about each of the following situations:

1. A fixes a bug.
2. A adds functionality without disturbing existing functionality.
3. A doesn't respect the concept of "backward compatibility" and modifies EP in such a way that existing functionality is altered.

Your response to any of these situations could be (and I'm supposing that you always have access to the current EP source code through, for example, GitHub)

• A. to simply continue with the existing EP,
• B. to use A's new version of EP,
• C. to modify your copy of EP.

(There may be other options in each set that I forgot.) The Cartesian product of those two sets has 3x3 = 9 situations for you to consider. My point in this post is not to suggest how you should handle these situations; at this point, I'll leave that up to you. But I am saying that you should, at this point, decide how you want to handle each of those situations, because your decisions will have a major impact on the answer to the Question that you originally posed.

[These are just my own ideas that I've developed over the decades. Does anyone here know of a book or other source material where these ideas are discussed formally, something like Software Package Management for Dummies?]

It would be much easier and faster to recreate the data each time. You just want to ensure that you get the same random numbers each time. You can do that with option seed:

GenerateGaussian(100, 2, .6, seed= 42);

This will give you the same 100 numbers when run again, even if it's run on a different computer. If you want each student to have a different set of 100 numbers, just have them change the 42 to some positive integer of their choosing. When they submit their workbooks, and you re-execute them, you'll get the same numbers that they had. This option completely avoids the need for external files.

You haven't provided a definition for the function f. Thus f(X[i]) remains a symbolic expression rather than becoming an explicit number. Thus printf objects to you trying to use it in a position where it expects a number.

The error message refers to fprintf rather than printf, which could be slightly confusing. This is because
printf(...)
is identical to
fprintf(terminal, ...),
where terminal is simply the symbolic file identifier of your display screen.

The second derivative of f at 0 is (D@@2)(f)(0). What you had, D^2, is the square of the first derivative. The second derivative at 0 could also be expressed as (D(D(f))(0) or D[1,1](f)(0). This latter indexed form indicates specifically the second derivative with respect to the first variable and is the form that must be used when there are more than one independent variables.

Since 0.005 * 0.07 * 0.15 * 5000. = 0.2625 < 1, it is clear that the population is doomed. A discrete model of it is provided by:

#Returns fraction surviving after t time periods:
P:= proc(t, factors::list(positive)) 
local r;
   if t::nonnegative then 
      mul(factors)^iquo(trunc(t),nops(factors),r)*mul(factors[1..r])
   else
      'procname'(args)
   fi
end proc:

plot(P(t, [.005, .07, .15, 5000.]), t= 0..13);

And since we've multiplied the factors together already, it's obvious that a continuous model is simply

P:= t-> (.2625)^(t/4)

(Of course it's not worth it for this trivial example, but the efficiency of the discrete model could be improved by precomputing the partial products of factors. It'd be easy to have this done in a totally automatic and general way.)

If we let PP(p) denote the number of prime-inverse pairs in Z_p (using your definition from above), then there's strong theoretical reason to expect that PP(p) ~ Li(p)^2/(2*p). (I didn't look that up; I just derived it with a minute of thought.^[*1]) The computational evidence shown below supports this. Furthermore, the residuals seem to be bounded by (2/Pi^2)*sqrt(p). The sqrt(p) seems strong (from the computational evidence alone), and the 2/Pi^2 is just my guess. Here's my code and plot:

#This code uses Maple 2018. If you're using an earlier Maple, I'll gladly retrofit it. The only thing that'd 
#need to be altered are the 3 embedded assignments.

CountPairs:= proc(p::prime)
option remember;
local b, a:= 1, n:= 0;
   while (a:= nextprime(a)) < p do 
      if a <= (b:= modp(1/a, p)) and isprime(b) then n:= n+1 fi
   od;
   n
end proc:

#Compute as many as we can within a certain time limit:
try
   timelimit(
      999, #seconds
      #This infinite loop is intentional:
      proc() local p:= 1; do CountPairs((p:= nextprime(p))) od end proc() 
   )
catch:
end try:

#Extract remember table and convert to a Matrix:
PP:= Matrix(
   [lhs,rhs]~([entries(op(4, eval(CountPairs)), pairs, indexorder)]),
    datatype= float[8]
):
N:= max(PP):
Asy:= p-> Li(p)^2/2/p: #1st asymptotic term
Res:= p-> 2*sqrt(p)/Pi^2: #estimated residual bound
plot(
   [PP, Asy(n), Asy(n)+Res(n), Asy(n)-Res(n)], n= 2..N, 
   style= [point,line$3], color= [grey,red,blue$2], thickness= [0,1,0,0], 
   symbolsize= 1, symbol= point
);

As you can see, there are a handful of straggling outliers that are close to but not totally within the residual bounds.

^[*1]If we take as given that the events a::prime and (1/a mod p)::prime are independent, then the rest of the proof for the first asymptotic term is easy, and is left for the interested reader. The independence seems intuitively very likely, but I suspect that it's difficult to prove.

If you're willing to accept the handful of question-display formats that the package uses, then everything that you want can be done in a single line of code:

Grading:-Quiz("Remove the singularity from", x+1, (x+1)%*(x-2)%/(x-2), inertform);

Note the % signs used with the arithmetic operators. These prevent automatic simplification. If you're entering this in 2D input, you may need to use these inert operators in functional form rather than infix form. By giving the option inertform, this inert/active form distinction is 100% transparent to the student answering the question: They enter their expressions exactly as they'd enter any other Maple expression, but there'll be no automatic simplification.

If you're not willing to accept the question-display formats of Grading:-Quiz, then what you want (the prevention of automatic simplication) can still be handled by the InertForm package. 

You need to expand the numerator to prevent automatic simplification from cancelling the (x-2). For example,

Student:-Precalculus:-RationalFunctionPlot(expand((x+1)*(x-2))/(x-2));

Here's an example:

#List of desired numeric values of tickmarks (they don't need to be exact):
NTM:= [1e-5]:
 
plot(
   0, view= [(0..2e-5) $ 2], #just dummies for this example
   axis[2]= [ #axis[2] is the vertical axis
      tickmarks= [
         seq(
            #format is "actual numeric value" = "what you want printed"
            y= typeset(`1O`^nprintf("%d", ilog10(y))),
            y= NTM
         )
      ]
   ]
);

Details: ilog10(y) returns, as an exact integer value, the integer part of the base-10 logarithm of abs(y); in this case, that's -5. Then I used nprintf to to turn -5 into the string^[*1] `-5`. This hides that minus sign inside the quotes so that the two-dimensional processor typeset won't interpret it as something to put in the denominator.

By varying the above, almost anything that can be prettyprinted in a Maple worksheet can be similarly prettyprinted as a tickmark. You can even change the font, size, etc. Note that nprintf can take pretty much anything and turn it into a one-dimensional string of characters (even Chinese characters, for example, if you know their Unicodes); typeset handles the two-dimensional aspects, and is much more fussy and limited. More-complicated 2D cases require using the mostly undocumented Typesetting package.

Cases where you don't have foreknowledge of the precise numeric values of the y-axis tickmarks are much more difficult to handle. I've thought for years about that, but haven't come up with a general solution.

^[*1]In Maple-lingo, what I'm calling a "string" is properly called a "symbol".

The amazing thing about Maple is that in one short line of code, I can write a meta-procedure that broadly generalizes your assignment; and in a second short line, I can invoke that procedure to produce the procedure that directly fulfills your assignment:

restart:
MetaOp1:= (gen, n::nonnegint, pred, EA)-> ()-> EA(select(pred, ['gen()'$n])):
MyOp1:= MetaOp1(rand(1..100), 20, x-> x::even, x-> (S-> (S, S/nops(x)))(add(x))):

I'll show a test invocation of MyOp1 in a moment. Here, gen is a procedure that generates the numbers, n is the number of numbers to generate, pred is a predicate (such as Is x even?), and EA is the "external accountant" (by which I assume that your teacher  means the procedure that does the sum and average).

That doesn't use a for-do loop, so here's a meta-procedure that does:

MetaOp2:= proc(gen, n::nonnegint, pred, EA)
   proc()
   local R:= table(), k, r;
      for k to n do if pred((r:= gen())) then R[k]:= r fi od;
      EA([entries(R, 'nolist')])
   end proc
end proc: 
MyOp2:= MetaOp2(rand(1..100), 20, x-> x::even, x-> (S-> (S, S/nops(x)))(add(x))):

The arguments and invocation of MetaOp2 are exactly the same as those of MetaOp1. Now, the usage examples: In the following, randomize is used the get the same random sequence for both tests, but I've done it in such a way that you'll get a different random sequence each time that you run the pair.

randomize((r:= randomize())):
MyOp1();
                                 604
                            604, ---
                                 15 
randomize(r):
MyOp2();
                                 604
                            604, ---
                                 15 

I hope that you'll take this as instruction showing the broad possibilities of Maple and come up with your own solution to your homework.

It's very important for any Maple programmer to learn the trick (shown in MetaOp2) of accumulating results in a table within a do loop and then returning the entries or indices or both of that table. I'd probably fail a student who hadn't learned that by the end of the course. It's not always necessary to explicitly declare a table; for example, the L and L1 in Kitonum's Answer are tables. Never accumulate partial results in a list, set, or sequence; that's extremely inefficient.

I find this solution method to be much simpler than Kitonum's:

restart;
A:= <1, 4; 5, 1>;
x0:= <2, 1>; #initial condition
M:= Matrix(DETools:-matrixDE(A, t)[1]); #convert to Matrix
x:= M.(eval(M, t=0)^(-1).x0);

In particular, it avoids the need of the variables C1 and C2.

Note that I included an extra set of parentheses in the final line. Since Matrix multiplication is associative, these are mathematically superfluous. However, by using them to change the order of operation, the result is substantially simplified.