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.

Use

subsindets(p, patlist(0, anything), ()-> ());

This will remove any list (regardless of length) whose first element is 0. Of course, this is done without disturbing any other part of the structure.

I think that I can anticipate your next several questions, and I need to be away from the Internet for a few days. So, the code below will give you something to think about for a few weeks at the very least (and that would be if you study Maple all day, every day). But this code does seem to me to be the direction that you're heading.

The procedure below takes an unevaluated expression that can be evaluated at different levels (which I call a "cascade of evaluations" or "evaluation chain" or "evaluiation tree") and saves it as a vector of strings that can be used to reconstruct the entire cascade, even in a different Maple session, even on a different computer.

This uses some of the deepest non-mathematical concepts of Maple. It is a testament to the power and beauty of Maple that this can be done in such a short procedure.
 

Download NameCascades.mw

Like this:

algcurves[plot_real_curve]((x^2+y^2+2*y)^2 - 4*(x^2+y^2), x, y);

This only works for polynomials in two variables, but for such curves, this command gives much better results than other plotting commands.

Like this (is the simplest^[1] way):

MaxAbs:= L-> max(abs~(L));
and its usage would be
MaxAbs(L1);

You can (and probably should) add some type protection, like this (just an example of many possibilities):

MaxAbs:= (L::list(complexcons))-> max(abs~(L));

^[1]In some simple cases, such as this one, it is possible to create functions without using the arrow by building them from existing functions with functional operators such as @ (composition) and ~ (elementwise). So,

MaxAbs:= max@abs~;

When this is possible, it is often (but not always) more efficient than using the arrow.