@ivanfanthony 

I updated the code, and put in the parameter values that you gave. Please get the new code from the code window in the Answer above.

I think that there's an error either in your equations or in the parameter values, for these reasons:

1. The plots show extreme exponential growth. At k = 100, the values are up to 10^122.
2. The equations effectively contain (mu+sigma) and (mu+gamma), but using the values that you gave for those parameters, we have mu << sigma and mu << gamma, so the value of mu effectively contributes nothing to those computations.

@ivanfanthony I just posted my code with my ridiculuous values for the constants. You just need to change them in my code. The 11 lines where you change the values are all in one place.

Ideally, yes, we'd like K to be infinity. But the reality is that Maple cannot symbolically solve these equations. Nonetheless, the code I gave can theoretically compute any finite number of values. If you use a very large K and the computation is taking too long, I can easily make it run much faster (while sacrificing a little bit of the readability of the code).

@ivanfanthony 

Okay, give me some sample parameter values, and I'll give you code to plot them. I'll set it up so that you can easily change the parameter values. So, I need numeric values for these: N, S(0), R(0), I(0), E(0), beta, delta, gamma, mu, sigma, and the highest value of k.

@ivanfanthony So does that mean that I can just replace delta(k) with N?

@ivanfanthony You say that delta(k) and delta are both constants, but are they the same constant? That is, are they equal?

@ivanfanthony de

The S(k+1) equation has delta(k) while the I(k+1) equation has simply delta, which is inconsistent. So, is delta a function of k or a constant?

You say that you "can do the I(k+1) and R(k+1)". It's not possible to solve them separately. You either solve all 4 or none.

You need another equation, of the form S(k+1) = ....

It's unclear to me whether I is a function, a constant, or the specific constant sqrt(-1). If it's a function, then you also need an equation of the form I(k+1) = ....

Regarding your title: These are called difference equations. Differential equations are something else, containing derivatives of unknown functions.

@JoshuaLeiter Note that when Explore is used with parameters that cause an axis's nominal length to change (as the A parameter does to the y-axis in this example), then it's usually necessary (that is, if you want the plots to make sense as the parameters vary) to modify the plot command so that that axis's nominal length is the same for all values of the parameters. That's why my plot has the view option. The old-style animation commands (such as plots:-animate and plots:-display(..., insequence)) don't require this nuance.

@tomleslie Thanks. Yes, you're right, and I'm about to correct the Answer.

@acer In the case of my procedure / appliable module that's being discussed, the ModuleApply reconstructs/redefines/reassigns an exported procedure every time that it (the ModuleApply) is called.¹ It does this so that the exported procedure's remember table will start from scratch each time that the ModuleApply is called.² So, the remember table is not causing any time reduction for multiple calls to the ModuleApply; it's only reducing the time for a single call. Thus there's no need to use forget to get a more-accurate measurement of the time for an iterated call to the ModuleApply.³

Footnotes:
[1] I realize now that I could've accomplished the same thing by using a forget inside the ModuleApply. The only benefit of that would be making the code easier to read.

[2] And that is done because it's the correct thing to do for the algorithm rather than because it makes timing measurement easier.

[3] That leaves the issue of measuring the time required for the garbage collection resulting from the final (and only the final) destruction of the remember table. Of course, this amount of time is independent of the number of iterations, and, in this case, it's very small.

For almost any time measurement in Maple (not just the case at hand), there's still the issue of measuring the time required for the overall final garbage collection (meaning all garbage, not just remember tables). This time is, in general, not small enough to be ignored. None of timing techniques discussed so far handle this (although Acer did discuss the issue itself at length a few Replies above this). I have some ideas about this that are worthy of a separate post (probably several separate posts).

@sand15 There was a major bug in the most-recent code that I posted---just a two-character difference between what I tested and what got copy-and-pasted into the post. I don't know how it happened. I corrected the code in the Reply. You must take the new code. The matrices produced by the old code were totally garbage.

@Carl Love This next entry is 2 to 3 times faster than the one immediately above, so, at this point, it is by far the speed champion among all the candidates tested. This one uses somewhat-sophisticated table and list operations and very little arithmetic.

NimMatrix:= proc(N::And(posint, Non(1)), K::posint)
local 
   i, j, 
   NK:= N^K,
   baseN:= table(
      [seq(0..NK-1)] 
      =~ [seq([seq(i)], i= Iterator:-MixedRadixTuples([N$K], 'compile'= false))]
   ),
   unbase:= table((rhs=lhs)~([entries(baseN, 'pairs')]))
;
   `mod`:= modp;             
    rtable(
      'symmetric', (1..NK)$2, 
      [seq([seq(unbase[baseN[i]+baseN[j] mod N], i= 0..j)], j= 0..NK-1)],
      'subtype'= Matrix
     )
end proc:

@Ritti Why do want Array output? And why do you want an Array so small, 11x11? I would make a 3-D plot, play around with some options until I was satisfied with it, then, maybe, extract an Array directly from my favorite plot. Try this for starters:

sol:-plot3d(x= 0..1, t= 0..1);

Your seq command was too far from acceptable syntax for me to figure out what you intended. What range of x do you want? What range of t?

No Maple-related question is too simple here. But next time put it in the Questions section rather than the Posts section. I've already moved this one.

@Ritti The command is pdsolve. See help page ?pdsolve,numeric.