A list of lists, all having the same number of elements, is called a listlist in Maple-lingo. It's a primitive form of matrix, and it's often a useful substitute for an official Matrix. If A is a listlist, then its transpose is

`[]`~(A[])

Thus, what you want to do can be done by

A:= [[0, 1, 2, 3, 4, 5], [0, 2, 4, 6, 8, 10]]:
writedata(test, `[]`~(A[]), integer);
0	0
1	2
2	4
3	6
4	8
5	10

I don't build web pages, so I've never done this, but I suspect that you could embed a Maple Player  , which is a freely redistributable app whose main purpose is allowing those who don't own Maple to view your worksheets.

Maple can solve nonlinear ODE BVPs both symbolically and numerically. It can solve linear ODE IVPs by Laplace transforms, but not nonlinear ones.

Please post an example BVP of the type that you work with and state whether you want a symbolic or numeric solution.

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