The easy way to make a neatly formatted table is to use a Matrix:

interface(rtablesize= 21): #Allow display of matrices up to size 21x21.
Matrix([seq([t, evalf(-sin(t)), -.5*sin(.5*t)], t = 1 .. 20)]);

That's all there is to it! Note that I changed your curly braces {} to square brackets []. Don't use curly braces in any situation where the order of the contents within the braces matters.

With a tiny bit more effort, you can add column headers:

Matrix([[`time`, `Wave A`, `Wave B`], seq([t, evalf(-sin(t)), -.5*sin(.5*t)], t = 1 .. 20)]);

(which explains why I made the rtablesize 21 instead of 20).

The above can also be done with a DataFrame (see ?DataFrame), but, in this case, I didn't think it worthwhile to go that route.

The vector field <x, y> is not the same thing as the scalar function x*y. What you need is a path integral, not a line integral:

VectorCalculus:-PathInt(x*y, [x,y]= Path(<4*t, 3*t>, t= 0..1));   
     20

Or you may prepend Student:- to that; the syntax is identical.

Okay, here's a draft of a worksheet that you can start using right now. Make sure to heed the comment (set in a box of exclamation points) about reducing the pixels if your display is less than 4k.

I refer to your Q(x) as forcingfunc. I set it to simply x in this line:

OneFrame:= Explore_BVP(forcingfunc= x, scenes= [[[x,w(x)]], [[x,u(x)]]]):

You change that x to any function that you want. You can also arrange and combine the plots any way that you want (currently limited to horizontally) using that scenes option that I made up. There can be any number of plot windows left to right. Each sublist is one plot window. Each sub-sublist is a curve. Each algebraic element is a plot axis. Each window can be plotted in 2 or 3 dimensions depending on the number of elements in the sub-sublist. So, for example, if you want w and w' on the left plot, u vs. w on the middle plot, and a phase portrait for u on the right plot, you'd do

scenes= [[[x,w(x)], [x, diff(w(x),x)]], [[w(x), u(x)]], [[u(x), diff(u(x),x)]]]

If the refresh rate of plots is too slow as you move the slider, or maybe if you're getting a lot of popup error messages, change numpoints= 10000 to a lower number.

Download Explore_BVP.mw

(This one is a FAQ (not that we have a list of such).)

The keyword itself (i.e., the part that's on the left side of = when it's passed) is an unassigned global name. Thus, it needs an unary :- prefix to disambiguate it from the parameter of the same name (as is always the case when a global and some other name are ostensibly the same: Prefix the global). So, you need to do:

procB :=proc(the_equation)
  print("inside procB, the equation is ", the_equation);
  procA(':-the_equation'=the_equation); 
end proc:

The unevaluation quotes are optional (although good safe coding practice), but the :- prefix is mandatory in ambiguous cases. Some people use the prefix all the time; I find that excessive and ugly. If you don't already know that you've used an ambiguous name in your own procedure, well, then you're in pretty bad shape.

Note that no number of unevaluation quotes can ever turn a parameter back into its parameter name, as this example shows:

proc(x) print('''x''') end proc(5);

This is because the substitution of passed arguments for parameters happens before the execution of any code in the procedure's body.

You need to put y=5 in parentheses, (y=5), to disambiguate.

None of Maple's relational operators =, <>, <, <=, subset, in, etc., are associative: They all require exactly two operands; no more, no less; one on the left, one on the right. So, the expression A = B = C is ambiguous: Does it mean (A=B) = C or A = (B=C)? It's up to you to specify.

You may argue: "Of course = is associative, as well as symmetric and transitive!" Sure, that's the way that it's generally used in mathematical writing. But in computer languages, including your beloved Mathematica, it tends to also be used antisymmetrically, for example, as a means of passing arguments. Even in mathematical writing, "Let x = a" means something different from "Let a = x".

You've accidentally discovered that unevaluation quotes can act like parentheses. Try not to rely on this (unless you're code golfing). The quotes are very ephemeral. Examine the resulting expressions with lprint, and you'll see the parentheses.

What is Q? Even if you know what Q is, it's clear from your output shown above that you haven't told Maple what Q is. Or were you hoping to get a solution for an arbitrary forcing function Q? [LMFAO] Maybe that's what Maple is trying to do....

I suspect that if a symbolic solution is possible at all, it'll be enormously complex with deeply nested integrals and RootOf expressions using your C__s as coefficients. And what's the point of having an enormously complex symbolic solution that's incomprehensible to the human mind? I'll admit that there is sometimes a point to that, but what's your reason? Chances are good that a parametric numeric solver, specific to this BVP system, would suit your purpose better than a symbolic solution. That numeric solver could be written in 10-15 minutes, including on-screen sliders to change the solution's plot as the parameters are varied.

The numerical solution of the original ODE system can be obtained easily, quickly, and straightforwardly:

restart:
F:= [S,E,J,R]:
IVP:= [
   diff(S(t),t) = mu*N - (beta*J(t)+mu)*S(t),
   diff(E(t),t) = beta*S(t)*J(t)-(mu+sigma)*E(t),
   diff(J(t),t) = sigma*E(t)-(mu+gamma1)*J(t),
   diff(R(t),t) = gamma1*J(t)-mu*R(t),
   S(0)=S0, E(0)=E0, J(0)=J0, R(0)=R0
]:
params:= [
   N= 5e7, beta=0.00001, mu=0.02, sigma=45.6, gamma1=73,
   S0=12500000, E0=50000, J0=30000, R0=37420000
]:
sol:= dsolve(IVP, numeric, parameters= lhs~(params), maxfun= -1):
sol(parameters= params):
plots:-display(
   <seq(
       plots:-odeplot(
          sol, [t,f(t)], t= 0..70, numpoints= 5000
       ),
       f= F
   )>, thickness= 3
);

But I don't know how to solve the system after you've applied "homotopy...".

In the code below, I used gamma1 for gamma and J for I because the symbols gamma and I are reserved in Maple for specific constants. If you really want to use gamma and I, it's easy to make a few adjustments to accomodate that.

restart:
Sys:= [
   S(k+1) = (mu*N - beta*Sum(S(l)*J(k-l), l= 0..k) - mu*S(k))/(k+1),
   E(k+1) = (beta*Sum(S(l)*J(k-l), l= 0..k) - (mu+sigma)*E(k))/(k+1),
   J(k+1) = (sigma*E(k) - (mu+gamma1)*J(k))/(k+1),
   R(k+1) = (gamma1*J(k) - mu*R(k))/(k+1)
]:
V:= [S,E,J,R]: #The unknown functions
ICs:= V(0) =~ [S0,E0,J0,R0]: #The symbolic initial values

K:= 100:        #highest value of k to compute and plot
params:= [
   N = 5e7,     #population
   mu = 0.02,     #?
   beta= 1e-5,     #?
   gamma1= 73, #?
   sigma= 45.6,  #?
   S0= 1.25e7,     #S(0)
   E0= 5e4,     #E(0)
   J0= 3e4,     #J(0)
   R0= 3.742e7     #R(0)
]:
assign(eval(ICs, params)); 
NumSys:= subs(Sum= add, unapply(eval(Sys, params), k)):
for _k from 0 to K-1 do assign(NumSys(_k)) od:
plot([seq([seq([k,f(k)], k= 0..K-1)], f= V)], legend= V, labels= [k, ``]);

Use this:

SP:= (n::posint)-> add(mul~(ifactors(n)[2])):
SP(12);

The [2] part remains the same regardless of n. It's an index into the list structure returned by ifactors. See ?ifactors, and note the (slight) distinction between it and the more commonly used ifactor. If you use ifactors  instead of ifactor, the format is more consistent, and you're less likely to encounter exceptional cases. The ifactor form is intended primarily for display, while the ifactors form is intended primarily for subsequent computation.

Your line

CC(i, j) := CC(i, j)+Testna(i, 1+(n-1)*3)*(Z(1, n)-Z(1, n+1))

should be

CC(i, j) := CC(i, j)+Testna(i, j+(n-1)*3)*(Z(1, n)-Z(1, n+1)),

the only difference being that I changed 1+... to j+....

Your whole thing could be made much simpler, and made to work for arbitrary K, like this:

K:= 4:
Z:= Vector[row](K+1, n-> z-2*(n-1));
DZ:= Z[..-2] - Z[2..]; # minus Delta Z
Testna:= LinearAlgebra:-RandomMatrix(K-1, K*(K-1)):
CC:= Matrix((K-1)$2, (i,j)-> add(Testna[i, j+(K-1)*(n-1)]*DZ[n], n= 1..K)); 

[The third line of the above code has been corrected due to an error discovered by Tom, below.]

So, there's actually no explicit loops required because the looping is built-in to the Matrix and add commands, and it's usually more efficient as well as much more convenient to rely on these built-in features.

The preferred form of Vector and Matrix indexing uses square brackets rather than parentheses. So, it's Z[n] and Testna[i,j] rather than your Z(1,n) and Testna(i, j). Your parentheses also work, in this case, but I think that it'd be safer for a beginner to stick to square brackets, as the parentheses form has a special meaning that I'd rather not go into here.

Using Maple's Standard GUI (which you're probably already using anyway), try this:

(xrange,yrange):= (-4*Pi..4*Pi, -2..2):
Explore(plot(A*sin(B*x), x= xrange, view= [xrange,yrange]), A= yrange, B= -2..2);

and play with the A and B sliders that'll appear on the screen beneath the plot.

This specific command doesn't animate, but I wanted you to "explore" it first before we animate it. And I think that any pedagogic purpose served by an animation is better served by the above command, assuming that your purpose is pedagogic. Explore has a vast number of options, including many for animation. See ?Explore.

You can use Diff as an inert form instead of diff. Replacing diff with %diff should've worked though.

When you call it an "operator", I think that you mean that you intend to "operate" on an unknown function f, making a new function. Is that correct?

It is done by using the output option of Usage, as in

Time:= CodeTools:-Usage(..., output= cputime);

A global variable and a local variable with the same apparent name are still completely different variables to Maple. Surely you know that already (at least on some level) because it's the same in virtually all computer languages. What's your reason for making C1 local to foo? I wouldn't try doing pattern matching with local variables at all. I could probably get something to work, but I think that it'd get rather convoluted.

The command that you need is subs (short for substitute).

eq:= x(t) = _C1*sin(sqrt(k)*t/sqrt(m)); #You don't need to type this.

subs([ _C1= A,  sqrt(k)= W*sqrt(m)], eq);

Note that the first substitution is rather obvious, but the second is more subtle. It has to be that way because the sqrt(k) and sqrt(m) are split in the original. The left sides of substitutions are rather limited (they need to be vaguely "atomic" ), but the right sides can be anything.

To other readers: Yeah, I know that I'm not using "atomic" in the Maple sense of that word defined at ?type,atomic.