@jgray07m You can read the help page ?Explore.

Can you come up with some numeric values?

@jgray07m Yes, it totally makes sense now.

There is really only one model because the first one is just a specialization of the second one at v=0. I guess that that's the vaccination rate.

The equations involving N are independent of the first three equations. Thus they can't change anything. So, why are they there?

To get any plots, we need numeric values, or reasonable ranges of values, for all 5 parameters, and we need 3 initial conditions (or reasonable ranges). I can make an Explore app where you can control the parameter values with sliders or knobs, and the plots will change with those.

@ The problem is on page 303 of the aforementioned book (by Steven Strogatz). I understand the differential equations now, and I'll enter them into a Maple worksheet.

@Rouben Rostamian  Nice work; Vote Up. I was on the verge of creating a similar graphic myself before I saw yours.

What you have presented inside the solve command is unintelligible to me as having any relation to differential equations. Where are the derivatives? Do you have two independent variables, phi and psi? There should probably only be one. Or is it an autonomous system with phi and psi as the dependent variables?

Maybe you should just try typing your equations here rather than in Maple. The copy and paste from Maple really mangled your expression.

There are no (ordinary) differential equations that I see in your worksheet. What are the dependent (or functional) variables? Are they H, S, R? What is the independent variable (there can only be one)? It's usually time t. 

You have 7 symbolic parameters (the lowercase letters). We'll need numeric values (or at least reasonable ranges of numeric values) of those in order to get a plot. We'll also need initial or boundary conditions such as H(0), R(0), S(0).

I see that you've computed a Jacobian matrix​​​​​​. That's not a usual step in solving a system of differential equations. 

Are you interested in doing it over finite fields of non-prime order? I could easily (an hour or so) write an object for you that does that. Can you restrict the field characteristic p to be greater than 3? That would be helpful because 2 and 3 need to be treated as special cases in elliptic curve arithmetic (and I'd need to do a bit more research on how to do the addition).

@HS In order for the "modding out by q" to happen in the background, you need to use any one of the 45 or so special commands that work with mod. These are listed on help page ?mod. One of those commands is Eval (capital E). So, if you do

Eval(f, x= q^2) mod p;

then the result will be reduced by q. It might seem from your example that simplify is "mod-aware"; however, it was just using q's definition as an algebraic number (from the alias) to reduce it over the field of rational numbers and then the integer coefficients got reduced by p. The command Expand (capital E) is one of those special mod commands, and can replace your use of simplify. In your if-statement, both sides of the equation can be put inside one Expand command:

if Expand(q^6 + q^2 + 1 = 2*q^2 + 2*q + 1) mod p then ...

You can easily add new commands to those 45, which is exactly what I did with Generators. I'll give the details of that when you want.

You never told me what Maple version you're using. I know it's not Maple 2019, so I know to avoid the newest syntax. But perhaps I need to avoid Maple 2018 syntax also.
   

@tomleslie The problem with this is that there are two real roots very close to x=0 for any a in (0,1] (degenerating to a multiplicity-2 root at x=0 for a=0). It's easy to prove this with elementary calculus or with solve(..., x, real, parametric) (although this latter technique is far from elementary "under the hood"). The original plots as well as yours decay near a=0, although yours get closer to 0 due to increased precision.

It's possible to get perfect real-root plots from the parametric solutions. It's just tedious to deconstruct the nested-lists-of-Rootofs-inside-piecewise form of the solutions. But I have done it in a somewhat automatic way, and I'll post it when I get to a computer. 

I don't know what you mean by "texture" fill, but I recall the fill algorithms available circa 1980 (Pascal or otherwise) being limited to convex regions. 

Please don't post Questions in inappropriate places. This Question has nothing to do with C-means clustering, so I've reposted it as a separate Question.

@jgray07m You need the proper Maple syntax for the fundamental constants. It should be

z1:= exp(Pi*I/3);  z2:= exp(-Pi*I/3);

So note that Pi and I are capitalized, and e is replaced by the exponential function, exp.

@HS I used syntax new to Maple 2019. What Maple version do you have? The following will work as long as your Maple has index (Maple 2017?):

`mod/Generators`:= proc(q::RootOf, p::prime)
    local 
        x, n:= p^degree(op(q))-1, 
        Ns:= n/~NumberTheory:-PrimeFactors(n):
    select(
        e-> not 1 in Expand~(e^~Ns) mod p,
        index~(Roots(x^n-1, q) mod p, 1)
    )
end proc:

Check if your Maple has index. If not, I'll need to modify the above.

The command to look up in help is index, not index~. The ~ is the elementwise operator, which you can read about on help page ?elementwise. Basically, it makes any command act on the objects in a container individually rather than acting on the container itself.

If your Maple version is earlier than Maple 2018, then I'll also need to change the until clause that appears later in the worksheet.

Regarding syntax highlighting, the Code Edit Regions in Maple 2019 have it. Maybe Maple 2018 also.

@HS I think that there's a possibility that you may be confusing finite fields of non-prime order (which are of course constructed with polynomials) with the rings of polynomials defined over those fields. GF is an efficient choice (although the user interface is clunky) for computations in the fields themselves, and mod works for this also--still efficient, but not quite as efficient as GF^[*1]. But GF (as far as I know) doesn't work with the rings of polynomials defined over those fields; mod does.

[*1] because GF avoids some unnecessary conversions from/to external format to/from the much-more efficient zppoly format. 

@HS You wrote:

• Your way  for computing Gcd over non prime order fields does not work since your not moding out by an irreducible polynomial of degree 2 it won't even recognize when your polynomial is 0.

What you wrote there is totally wrong. You don't know what you're talking about, and I suspect that you have an attitude problem that's blocking your ability to learn. Do you think that I faked the output in that Answer? Do you see any q coefficients of degree 2 or higher in my gcd? The "modding out" that you refer to happens in the background, so you don't notice it.

Here's an example of computing generators and inverses in GF(7,3):
 

Download GF73.mw