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

@Christian Wolinski So, you managed to put the OP's code into an executable format, but you didn't think to upload it so that others could continue the work without needing to start from scratch? Collaborate.

@asma khan If matrices A and B have the same number of rows, then they can be attached side-by-side as Tom showed:

<A|B>

If they have the same number of columns, then they can be attached top-to-bottom by

<A,B>

These commands are not limited to two matrices. They can be interchangeably nested to arbitrary depth to create matrices of arbitrary complexity.

@nm I do not "object" to you saying "internal" error. I merely thought that you'd misinterpretted the int in `latex/int` to stand for "internal". On the other hand, with your clarification on your definition of "internal error", it seems to mean the same thing as "bug".

As usual, I'm just trying to offer a practical solution so that a user can proceed with their work. I'm not trying to excuse the error. I'm not addressing what should be; I'm addressing what is. I do think, like you, that the command "should  internally handle these issues", but I don't work for Maplesoft! I don't have the power to fix bugs! You should file an SCR. You and I have had this conversation several times before. Please do not respond to my workarounds as if I had the power to make things work as they "should". It is very annoying to me when you do that. I am just a "user" who likes answering questions and providing help.

Practical means, among other things, that the solution may not work on 10,000 equations all at once without error. I've never heard of anyone performing a large project with the latex command without needing to do some "hand tuning", and I'm amazed that you've not had a problem with the command before. It's a very weak feature of Maple, IMO.

If you'd like to correct the bug before someone at Maplesoft gets to it, the procedure `latex/int` would be trivial to correct by taking out the first two lines, which check that the integral has exactly 2 arguments.

@David Sycamore Got it. It wasn't clear from your original whether primality or parity should be checked first (and it only makes a difference for n=2). Also, you wrote a(n+1) = n+2 where it should be a(n+1) = a(n)+2. So with that in mind, the procedure can be simplified as follows. I also show one way you can get the pairwise output that you want. There are many ways that this could be formatted.

restart:
A:= proc(n::posint) 
option remember; 
local a:= thisproc(n-1)+1;   
   a +`if`(a::prime, 0, 2-irem(a,2))
end proc:
A(1):= 1: 
:
seq('A'(n)=A(n), n= 1..20);
   A(1) = 1, A(2) = 2, A(3) = 3, A(4) = 6, A(5) = 7, A(6) = 10, 
   A(7) = 11, A(8) = 14, A(9) = 16, A(10) = 17, A(11) = 20, 
   A(12) = 22, A(13) = 23, A(14) = 26, A(15) = 28, A(16) = 29, 
   A(17) = 32, A(18) = 34, A(19) = 36, A(20) = 37

And there's no reason that such a list must start at n=1. You could just as well do n= 100..120.

@mmcdara Yes, I'd say that it's never wise to use an indexed name such as x[0] together with its parent x when a subscripted name such as x__0 will do the job. The main problem arises when assignments are made. Any assignment made to x[...] when x is unassigned will turn x into a table. Since tables have the last name evaluation property, it'll still be possible to some extent to use x as a variable. Procedures also use last name evaluation, which is why the OP's use of Zeta doesn't lead to total nonsense. On the other hand, an unwary assignment to x will cause x[...] to be nonsense. 

I'm not on my computer right now, just on my phone, and I don't know much about such problems anyway. Yet your r < a sticks out to me. Shouldn't that be r <= 1, if anything at all? And how about r >= 0?

@mmcdara Your solution works for expanded polynomials (which is the case for this Question), but fails for almost anything more complicated such as rational functions or transcendental functions. Besides, there is a specific command to handle this precise situation gracefully: evalc.