What do you want it to return for x*y?

What does "2" mean in your (dx/dt)2?

@Carl Love I also attempted a quadratic residue analysis to prove that there were no other solutions. This, however, failed: I was easily able to find at least one quadratic residue (other than f(1), of course) for every modulus up to 2^17 = 131072. Here's my code:

F:= rsolve({f(n) = 4*f(n-1) - 2*f(n-2), f(1)=1, f(2)=28}, f(n), makeproc):
upperN:= 2^13: upperM:= 2^17:
FL:= F~([$2..upperN]):
for m from 2 to upperM do
   if andmap(x-> numtheory:-quadres(x,m) = -1, FL) then print(m); break fi
od:

@acer This is a great Answer, spectacular in several senses of that word!

A Post that I wrote (which you made significant contributions to) within the past year delves significantly into dropped-shadow contour plots (such as yours above) generated from parameterized  dsolve(..., numeric, method= bvp) output: "Numerically solving BVPs that have many parameters". (I know that you're already aware of this; I mention it for the benefit of other readers.)

@Jjjones98 The code above is a procedure definition; it is not intended to have any output. To get output, you need input, specifically the matrix size; you need to enter something like UMG(5).

@mmcdara 

In my opinion, those extra characters that may appear at the beginnings of continuation lines within an execution group that look the same as the prompt characters are just garbage, bugs, that only occur in some versions of Maple. They serve no useful purpose, are a total nuisance to remove, and I don't consider them to be prompts because a prompt is what your cursor sits to the right of while awaiting input. In other words, I'm saying that there's no such thing as a "passive" prompt; there are only active prompts and garbage characters.

However, I understand your point that my parenthetical comment could be misunderstood because those garbage characters do look similar to prompts; however, they can be distinguished from true prompts by paying attention to the fine, black bracket lines in the left margin.

@David Sycamore Kitonum's code checks the square-free property by checking that all exponents in the prime factorization are 1. The commands sqrfree and isqrfree are not for checking the condition. Rather, they are for producing a square-free factorization. This is a bit more elaborate, and it'd be wasteful to use these for the purpose of checking the condition.

@David Sycamore I am working on some code for this. I'll post in several hours.

It's far more efficient (and quite easy) to generate only the odd square-free integers than it is to generate the odds and then select those that are square free. This is why I didn't include a square-free check in my Answer.

Note that any odd square-free multiple of 3 will have the property that you're iinvestigating.

I deleted your separate followup Question, which seemed to be about this same code as above with some minimal corrections. Please post your followup Question here in this thread. I would've moved it here myself, but I don't have a way to do that and have it appear under your name.

Also, your Question is about Maple alone. It's obnoxious to check all the boxes: Maple, MapleNet, Maple toolboxes, MaplePrimes. Your Question isn't about any of those other things. But you should select a Maple version number.

Also, you should respond to Tom Leslie's Answer below. It's rude to not respond to people who've put significant effort into answering your Question.

Here's one possibility for implementing Dr Lopez's suggestion for your specific problem:

DynamicSystems:-DiscretePlot([seq(3*k+2, k= 1..6)], 1, 1, style= stem, view= [0..6, default])

There are many possible variations.

Would you please attach your worksheet as a downloadable file? Then we don't need to retype your equations. Use the green uparrow on the toolbar of this editor. Another acceptable format is to simply post your code in plaintext.

Just eyeballling it, it looks like you have a polynomial system that has total degree 36 in the decision variables and 12 symbolic coefficients. So, if any solution can be obtained, it's likely to be extremely complicated and abstract. It'll likely be in RootOf form (abstract algebraic numbers).

Edit: I see that Acer has solved it, and the solutions are not inordinately complicated. The estimate in the previous paragraph was just a quick mental computation, not accounting for unexpected simplifications specific to this system of equations.

@nm The workaround that I originally proposed

latex(subs(AllSolutions= (), sol))

works just as well whether sol is a single equation or a list of 10,000 equations. There'd be no need  "to check each time for something and make a workaround to make latex work."

@Christian Wolinski I used my Moderator powers to attach the executable--and now executed--code to the original Question. I also removed it from your Answer because you may not want to be associated with it, which is understandable. But feel free to put it back in your Answer if you want or to tell me to do so. Thank you for your cooperation and effort.

I've given you a vote up just for the painstaking work that you put in to unravel that "brick".

@tomleslie You're right, Tom, and I've now edited the executable code into the original Question.

@Christian Wolinski Thanks, and I understand your concern about not wanting to be associated with a "heap". (I tend to call it a "brick" of code :-).) Any executable form is better than the original. Sorry if I sounded harsh.