The numbers you quote for Mathematica (ie 10^308/10^-308) represent the standard range for double precision numbers, which is what I'd expect to see when using hardware floats. The numbers I posted above are the limits for software floats in Maple

I don't know anything about Mathematica, but doesn't it have a 'software-float' capability?

@Allan Wittkopf 

That didn't happen for me on 64-bit Linux in Maple 2016, but I don't have windows, so I cannot check. Also I'm not 100% sure I have the correct system, as I had to copy/paste the system from earlier in these threads, then manually edit for line breaks, as there were things like:

And I no longer have a dual-boot, so I guess we can't do a like-for-like comparison:-(

You can pick up my (functional?) worksheet diffEqProb2.mw, from one of my earlier posts on this thread. This is syntactically correct, and the only change you will have to make is for the 'path' to the MC.m file which actually contains the OP's equations.

Win64 is a bit of a special case for compile - you need to install the Microsoft visual studio compiler and the Windows SDK to have the compile option work there.

I have not used the 'compile' option for this problem: partly because the OP wasn't using it (so I didn't want to introduce another variation), and partly because the whole worksheet runs in more-or-less acceptable time uncompiled - so it didn't seem worth it. I have just added some timing points to the worksheet: the whole section utilising the lsode method runs in 4.7seconds, and the section running the rkf45 method runs in 65secs.

As for the maxfun, just set it to 0. This *should* work for you though I do not have Win64 to check:

dsolve(MC, numeric, method=rosenbrock, maxfun=0);

If you have been reading this - you will realise that it is all about setting "appropriate" values for maxfun stepSize, abserr etc parameters: your original suggestion was

sol := dsolve(MC, numeric, method=rosenbrock);

so I assumed that this would work - but, as reported earlier - it errored on a maxfun limit. Now could I fix this: probably. But I didn't want to get into fiddling with the options on yet another ODE solver.

However in the spirit of fair play, I swapped in the command dsolve(MC, numeric, method=rosenbrock, maxfun=0). The good news is that it runs with no errors, and gets the correct results. The bad news is that it takes 4477 seconds - compare with 4.7 seconds for lsode[adamsfull], and 65 seconds for rkf45. All times are for running without the 'compile' option.

Note though that if you are running into the maxfun limit, this suggests that dsolve is working pretty hard to solve the problem - is there any chance that the solution is near a singularity, and the round-off difference between Linux and Windows is sufficient to make that large a change? Do you know how stable the solution is?

Well the problems in obtaining the solution which have been reported on ths thread are indicatins that this is deeply unpleasant ODE system. If you actually look at the system from the MC.m file, you will see that one of the ODEs is lengthy, with multiple piecewise components, so I'm guessing it probably has no zeroth-order discontinuities, but probably quite a lot of first-order ones

@Allan Wittkopf 

Whilst I keep admitting that my solution is crap - it works

I substituted your suggestion for the first dsolve() command in my posted solution.

On Maple 2016, 64-bit Win 7, it errored out with

Warning, cannot evaluate the solution further right of .63477096, maxfun limit exceeded (see ?dsolve,maxfun for details)

I'm getting tired of criticism on this: which part of the word "works" don't you understand?

@vitato

The last comment I posted on this subject was

"Before I do any thinking about this - are you absolutely sure these equations are correct?????"

According to what Preben has posted - seems like your answer should have been - NO.

So far as I can tell, one solution for this problem has been posted.

Mine.

Prior to that solution posting there was a lot of basicaly unhelpful chit-chat After that solution posting there are a lot of "suggestions" for improvement, detailed explanations, whatever: but strangely enough - no other solutions.

Now if you read the soluttion I posted, there is a comment which actually says

#
# A really ugly way to get the the list of values you require
# This one could get me barred from MaplePrimes forever!
#

so trust me: I know the solution I posted is pretty crap. However

1. I posted a complete (if ugly) solution
2. No-one else has posted an alternative, so
3. Give me a break, because I won't be revisiting this unless someone posts a suitable, complete, alternative, working across both solution methods and in both Maple 2015 and 2016
4. So far the OP has not claimed that my seriously ugly solution is non-functional

@asa12 

You have three first order equations in three unknowns.

You therefore need three initial/boundary conditions

You cannot find/derive initial/boundary conditions from the differential equations themselves. They are extra information about your system which you have to supply. See any textbook on differetnial equations.

@vitato 

Now I know the equations to use, thought I'd do a quick check using

sys := {((D@@2)(phi))(t) = -2*(D(phi))(t)*(D(theta))(t)*cos(theta(t))/sin(phi(t)),
        ((D@@2)(theta))(t) = (D(phi))(t)^2*cos(theta(t))*sin(theta(t))-9.8*sin(theta(t))
       };
ICs:={theta(0) = (1/2)*Pi, (D(theta))(0) = 0, phi(0) = (1/2)*Pi, (D(phi))(0) = 1};
sol:=dsolve( (sys2 union ICs, numeric));
plots[odeplot](sol, [t, sin(theta(t))*cos(phi(t))], t=0..10);
plots[odeplot](sol, [t, sin(theta(t))*sin(phi(t))], t=0..10);
plots[odeplot](sol, [t, cos(theta(t))], t=0..10);

which should plot x(t), y(t), z(t), from 0..10. However Maple's numeric ODE solver complained of a singularity at t=~0.608. Seem rather unlikey that the equations of motion for a spherical pendulum should contain a singularity. Before I do any thinking about this - are you absolutely sure these equations are correct?????

The two definitions of your system are different:

sys := {((D@@2)(phi))(t) = -2*(D(phi))(t)*(D(theta))(t)*cos(theta(t))/sin(phi(t)),
            ((D@@2)(theta))(t) = (D(phi))(t)^2*cos(theta(t))*sin(theta(t))-9.8*sin(theta(t))}

from the dsolve() command

            ((D@@2)(phi))(t) = -2*(D(phi))(t)*(D(theta))(t)*cos(theta(t))/sin(phi(t))
            ((D@@2)(theta))(t) = (D(phi))(t)^2*cos(theta(t))-9.8*sin(theta(t)),

Emboldened term in first definition does not exist in the second!!

Which is correct????

@asa12 
Think about the statement

dsolve({eq2a = exp(t), eq3a = exp(t)});

Since eq2a = exp(t) and eq3a = exp(t), then it must be true that eq2a = eq3a, which is a single equation with two dependent functions: so I assume that you accept that this cannot be solved explicitly. The best you could hope for is to get one of these functions in terms of the other. And since you now know why it is impossible to solve explicitly, you would not expect Maple to provide such a solution - would you?

This is essentially what the dsolve() process is doing. There is no possibility of explicit functions for b(t), c(t), so  the best it can do is to return relations between the two functions;

@adrian5213 

Christopher2222 is correct. I have amended the code in the attached worksheet

randBin.mw

1. You define eqn1, but then use equ1: fixed this, but then
2. Executing indets([equ1, equ2, equ3, equ4, equ5, equ7, equ8, equ9]); indicates that there are 13 unknowns and only 9 equations!!
   

"I found the command as ilaplace" -I assume you mean MTM[ilaplace](), which is actually the inverse transform. the forward laplace transform would be MTM[laplace]

Why are you using these from the MTM package?

The more conventional approach would be to use inttrans[laplace] and inttrans[invlaplace]

@marc sancandi 

"I guess maxfun=-1 in Maple 2016 is intended to be the equivalent of maxfun=0 in Maple 15 (?) : it authorizes  an unlimited number og function evaluations ... Is that it ?"

Must have had a brain-fade - I should have been using maxfun=0: however I have just tried this and it makes no difference. I still get the "excessive amount of work error". I have just been re-reading the dsolve help(), in particular the page

?Option maxfun for Numerical Solution of ODE Initial Value Problems

Now on this page it is claimed that "For the dverk78, lsode, and gear methods, this option is disabled by default." So it is already disabled! So "disabling" it some more isn't going to help. On the other hand the help page ?dsolve/numeric/lsode/advanced does list a parameter which controls the maximum number of function evaluations, with no indication of how to disable it.

Skip to the second par of your answer. I too observed that reversing the order of point evaluation (her back from 8.3s to 7.8s) changes the results ... and (I am not sure always) retrieve the ones corresponding to the odeplot

I thought I had explained this reasonable clearly - obviously not! I'll have another try. In order to get a "point" solution, at (say) t₁ then Maple has to run the solution procedure from the Initial Value Point to t₁. However for the next point t₂ (assuming abs(t₂) > abs(t₁) ), it only has to run from t₁ to t₂. This is one of the things which makes the odeplot() command fast.

However if you ask for a point solution at t₁ followed by t₂ where abs(t₂) < abs(t₁) then Maple has to compute this by running the solution procedure from 0->t₁, and then 0->t₂ Computing points "from "scratch" in this way can get very expensive and is more prone to stepSize related errors (which might show up as "different answers")

4/ back to your first claim (I apologize for this messy answer) "The fundamental problem is the default stepSize/error control which is being used." : I would me puzzled if the command "sol := dsolve(...)" ende in some error because of a bad time step or tolerance, or an insufficient amount of maxfun (see the test with rkf 45 where the default "maxfun=30000" has been set to 500.000 to obtain a full range solution [I did not check for sol('numfun') to see the real number of evaluations rkf45 did]).
I often face these kind of difficulties (which is the reason why I test different solvers before using one for a specific kind of problems)

It seems reasonable to think that if the dsolve procedure ends without error and reaches the highest time of interest (say here 10s) ; and the odeplot is (seems) correct ... than I expect a point evaluation is correct.

If you were making ths point for functions with continuous first (and preferably also second) deivatives, I might be tempted to agree with you. However if a DE solver manages to "step over" a discontinuity then all bets are off. It may generate a completely erroneous answer. I managed to achieve this in my previous post with the rkf45 method, where the output was fine up to t=~3, but then "flatlined" all the way to the end in the odeplot() command. I assume some "event" happened in your system at this time, and the rkf45 method on default accuracy missed it - then all subsequent data is garbage. Now I managed to fix this just by crunching 'abserr',
 but there are other ways to do this - see the ?Events for dsolve[numeric] help page. Utilising techniques on this page can greatly help all/most of the DE solvers, when there are known first (ir second) order dicontinuities in the problem

The fact that this can occur in odeplot() with no warnings/errors using the rkf45 method should help you appreciate that a DE solution procedure can occasionally return garbage results. Obviously this can also happen when seeking a "point" solution.

@Carl Love 

Interesting

?dsolve/numeric/lsode
?dsolve/numeric/lsode/advanced

Reading these two pages together it is fairly obvious that the underlying lsode method described on the second has a *different* interface from the other ODE solvers - references to 'atol', 'rtol', 'mxstep' etc rather than 'abserr', 'relerr', 'maxstep' - also the fact that these are to be entered as numerical values in a vector, and their purpose is interpreted according to position in that vector.

I'm guessing that the first of the above pages is describing a simple "wrapper" which makes the interface to this method *look* like other Maple solvers (same option names etc), but all that really happens is that this "wrapper" just constructs the appropriate control vector for the underlying method.

The maxstep/mxstep issue

?dsolve/numeric/lsode lists 'maxstep' as an option, and this is what I used (and it worked)
?dsolve/numeric/lsode/advanced lists 'mxstep', as an option.

However their purpose is different!

'maxstep' does stepsize control, and presumably gets mapped c₇=hmax

'mxstep' looks as if it is doing roughly what 'maxfun' does, and is presumably mapped from it

The description of the remaining error I was getting as "an excessive amount of work" exactly matches the description of an  'mxstep' problem.  - so maybe the 'mapping' from 'maxfun' to 'mxstep' doesn't work when I try to turn it off completely with maxfun=-1 ????

I might try playing with the underlying (more direct interface), but I doubt it! I'd have to do too much reading, to figure out values for some of the parameters esp c₃, c₄, c₅

I think what we all want is a worksheet we can "play" with.

The attached is my attempt to provide this: it comprises

1. The OP's original code with modified read statements to import the relevant differential equations
2. Note that you will have to download the MC.m file from the OP's specified dropbox location. Depending on how you handle downloads - you ought to know where this downloaded file will appear
3. You will then have to modify a couple of 'path' statements in the attached, to match your specified download directory. I have color-highlighted the relevant parts
4. With these mods, I now have a worksheet which actually executes, although I have performed no analysis/diagnosis. I'm pretty sure sure it replicates the OP's original issue, but have not competely checked
5. I have just poured myself a rather large scotch, and bedtime is approaching, so any further response may take a while

diffEqProb.mw