@EvilDuncan The answer changes for me. Using HeffR, I have the eigenvalues expressed as the roots of a relatively simple polynomial. That is, its coefficients are relatively simple compared to the algebraic functions that were input. You can also use specific values of x as in Eigenvalues(eval(HeffR, x= 1)).

Can you post an executed worksheet that shows them unchanged?

Ultimately, you're going to need decimal approximations for these eigenvalues, right? When you get to that stage, you'll get more-accurate results from

Eigenvalues(evalf(eval(HeffR, x= a)))

(where a is some explicit numeric constant) than you'd get from decimal approximation of the roots of the polynomial.

@Adam Ledger Ah, but there is a field of math called "catasthrope theory". I even took a one-credit course in it once. That was a very long time ago (40 years), and I don't remember much. Anyway, see the Wikipedia article "Catastrophe theory".

@Fabio92 Fabio, you're missing a colon in the assignment to pw.

@maple2015 Good point about the first-order/second-order. For numeric solution, dsolve converts all ODEs to first-order systems (in the usual most-trivial way).

Did you try the stepsize option? And you still got a too-accurate solution? The classical methods are not supposed to use "controller steps".

I need to go to sleep now. 5:14 AM here. I'll look at it more when I wake. I'll also explain how to dig through the code.

@inteli I found the problem. It wasn't your fault. In the original lines, the Arrays need to be converted to Matrices. So, add Matrix to all the Import lines:

M1A:= Matrix(Import(...));

etc.

I don't know why Tom didn't have the same problem.

@maple2015 By default, dsolve is using h=0.005 (according to the help page). That's much smaller than your h, which is why there's no visible error. You can control it by adding the option stepsize= h to the dsolve command.

In principle, yes, you can dig through the dsolve code and find the algorithm, but it's nasty work in the case of dsolve. The formula used for adambash is given on the help page as

Y[n+1] = Y[n]+(1/24)*h*(55*f(t[n], Y[n])-59*f(t[n-1], Y[n-1])+37*f(t[n-2], Y[n-2])-9*f(t[n-3], Y[n-3])).

Perhaps that's all that you meant by "the correct algorithm which is used by Maple".

@Preben Alsholm Thanks, I wasn't aware of that; I've always added square brackets. But it makes sense; many kernel commands have similar behavior---the ability the treat non-explicit sequences as single entities.

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