If your expression containing theta is expr, then do 

diff(subs(theta= theta(t), expr), t);

In other words, theta has been replaced by theta(t), which indicates that theta is a function of t.

Here is a simple closed-form answer that only uses real numbers:

U := n->
   exp(.238022513478555*n)
    * (1.70226106278932*sin(1.91811758317245*n)
       + 1.28091347285205*cos(1.91811758317245*n))
    + 6.71908652714795*exp(.622567261710995*n);

To get exact integer values, apply round.

I am working on an elementary expository derivation of the above formula for you. 

Your initial conditions are wrong. They should be U(-2) = 2, U(-1) =2, U(0) = 2.

What I present here is not the standard straightforward explanation; rather, I think it's the most-compact elementary explanation. Nonetheless, you should be able to find it in any good textbook for a full-year introductory single-variable calculus course.

We'll suppose that a(x) > 0 for all x. This is only done pedagogically to avoid imaginary numbers; it's not mathematically necessary.

Let y(x) = a(x)^b(x)  =>   ln(y(x)) = b(x)*ln(a(x)) . Differentiate both sides wrt x:

y'(x)/y(x) = b'(x)*ln(a(x)) + b(x)*a'(x)/a(x)  =>  y'(x) = y(x)*(b'(x)*ln(a(x)) + b(x)*a'(x)/a(x))  

=>  y'(x) = a(x)^b(x)*(b'(x)*ln(a(x)) + b(x)*a'(x)/a(x)). Now, rewrite "prime" notation as "diff" notation and you'll have what you posted.

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