@radaar To be even more specific topologically, do you mean embeddings^[*1] of [0,1]x[0,1] into R³? So, for example, a mountainous region of the Earth's surface, but not the whole Earth itself?

[*1] An embedding f:A->B is a bicontinuous injection, i.e., a homeomorphism between A and its image in B, (f(A)).

@Mariusz Iwaniuk 

The final condition is D(f)(infinity) = 1, not f(infinity) = 1.

When you put comments in worksheets, they should be marked as comments so that the syntax checker does not interpret them as code, which'll surely produce error messages. There are a few ways to do this. The easiest is to begin the comment with #. Then everything to the end of that line will be ignored by the syntax checker.

It's also possible to set aside whole execution groups (which are somewhat analogous to paragraphs) as fully formatted text, similar to that produced by a word processor. 

@Zeineb 

X_k is the k^th column of A^-1.

@jgray07m So, populations are relative, as in percentages, rather than actual population numbers? Yes, it seems reasonable. I guess then that N is total (alive) population? 

And I need the 5 parameters. Obviously, the v slider should go from 0 to 1.

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