@Kitonum Your display command is missing seq.

@Sabrina Kamal Use

plots:-display([seq(rect[i], i= 1..26)], color= black);

Some comments on your programming style:

1. A lengthy formula (say, one that's wider than 1/3 the screen) should never appear in your worksheet more than once.
2. You should never copy-and-paste the output of one command as the input of another.
3. You should compute the error for each value of x and plot it, perhaps on different axes. Also, compute the maximum absolute value of the errors (use max command).
4. You should do many values of h. Make a plot and/or a log plot and/or a log-log plot of max error vs. h. It should be a fairly smooth curve, likely a power function. Use Statistics:-PowerFit to estimate it. The exponent on the h is the important thing.

I can help you with any of these.

@Sabrina Kamal At the end of your most-recent worksheet, you simply need to add the command

plots:-display([seq(rect[i], i= 1..26)]);

@Christian Wolinski Why don't you experiment with it to find out? It's the limitations of your system that matter, not mine.

@Christian Wolinski fsolve can count the real roots very quickly, but there may be a very few real roots that it doesn't find. Set Digits to a high value. It'll eventually find them all if you use enough Digits.

RootFinding:-Isolate claims to be able to isolate real roots exactly (into intervals with rational endpoints). I don't have experience with it on a polynomial of this size to know how fast it is.

I think that there may be a way to embed the real symmetric matrix into a larger matrix, preserving its real roots, and such that the characterictic polynomial is unlikely to split.

@Christian Wolinski What dimension do you want to increase to find the limits?

Increase the 600? I think that you can increase this only a little.

Increase the 200? Yes, I think that you can easily increase this into the thousands.

Increase the {-1, 0, 1}? I think that it'll get bogged down from just a small increase.

Increase the density 1/200? You can increase a little before getting bogged down. I obtained best results using density= 3/n with n=600. If you lower it to 2/n, then 0 will appear as a root many more times. Of course, it's trivial to factor out any discard 0 roots.

Notice that the two ODEs are not coupled. If you separate them into two separate dsolve commands, you'll get solutions. For one of these, you'll need to apply value to force it to "do" the integration, even though it has no trouble doing this. You might as well just apply value to both. (Using value on something that doesn't need it simply returns the argument, so there's no harm.)

Here it is in a downloadable worksheet.
 

Download rectangles.mw

@inteli The error messages combined with Tom's worksheet shows me clearly that you're trying to append a column with 7 entries to a Matrix with 6 rows. Obviously that can't be done. Tom's worksheet shows that you have 2 data sets, one with 6 rows, labeled JAN through JUNE, and one with 7 rows, labeled JUNE through DEC. You must be inadvertently trying to combine one with the other.

Note that JUNE occurs twice. Is that a problem? The data for the two JUNEs are different.

Would you like to combine the two data sets into a single data set with 13 months? That can be easily done.

If you continue to have problems, please upload your worksheet that shows the error message. It only needs to be executed up to the first error message.

@nm I agree wholeheartedly; I hate reading superfluous parentheses in standard algebraic expressions.

Maple's conversion of 2-D input to 1-D form also suffers from automatically added superfluous parentheses.

Your Mathematica example doesn't look so great to me! The spacing is awkward, and the 2 is not superscripted.

@Sabrina Kamal It's not some unprintable character. The function's name literally has a question mark. Explicitly, it's

back-quote question-mark left-square-bracket right-square-bracket back-quote

Or, if execute this, you'll see the name in Maple:

convert(cat(StringTools:-Char~([96, 63, 91, 93, 96])[]), name);

@Sabrina Kamal In that case (Maple 2016) you can use

R:= unapply(map(`?[]`, <P1, P2, P3, P4, P1 | Q1, Q2, Q3, Q4, Q1>, [i]), i):
plot([seq(R(i), i= 1..21)], color= plots:-setcolors()[1]);

And before you ask, this neither goes inside nor outside your loop. Instead, it completely replaces your loop.

@one man Your solution works fine when the numerator of the fractional exponent is an even integer, but not if it's an odd integer. Your solution is probably fine for the OP in this case, but a surd solution is more general.

This doesn't seem to me like a problem appropriate for Groebner.

Do you want all solutions? Or would you be satisfied with just one solution? Are all the inequalities linear?