@momoklala 

I'm not psychic. You'd have to show me what the sytem looks like now.

When you report a problem, try to get into the habit of uploading a new worksheet. In this case it would be helpful if you just noted in your request for help, what change you have made. That would save me reading/comparing two worksheets to find out what is different this time!

@Preben Alsholm 

Have to agree :-(

In my own defence I was so pleased to get the rkf45 method working that I didn't check the solutions too carefully - the two methods agreed up to the point where the rkf45_dae method stopped working, and that was good enough for me!!!

I can only assume that the step size algorithm in the rkf45 managed not to hit 2 exactly. When a numeric solver "steps across" an undefined point such as this, then subsequent results will almost certainly be wrong

I note that the undefined point requires not only that U(x)=0, but also that x=2*x0. It may be generally true in the symbolic solution that U(x)=0, for x=2*x0, whatever value is selected for x0. I can't tell from the form of the symbolic solution which you supplied.

I guess I'll have to settle for the fact that the original ODE system can be transformed in such a way that the OP can now write his own rkf45 method, which seemed to be his objective. (Even if it will fail at x=2, if x0=1!)

@patient 

Before my last post I did a little digging into the diference between rk45 and rkf45_dae and couldn't find the distinction!!

Which is why I said "I do not know what the extension of the rkf45 method to rkf45_dae actually is"

Googling rkf_dae just gave mainly references to Maple! which we know does not explain the difference. This is why I made the comment "I would have to do a lot of reading to find out. (On the other hand you could do the same reading.)"

Just because you want a quick simple answer to your question, does not mean that one exists :-(

On the other hand if you type "dae numerical methods" into your favourite search engine you will get a lot of reading material

@patient 

No - the rkf45 method has not solved your problem!!!!!

What solved your problem was

method=rkf45_dae

If you try it with

method=rkf45

then you will receive the error message

Error, (in DEtools/convertsys) unable to convert to an explicit first-order system

If you check the help page for

dsolve/numeric/DAE_extension

then you will see the comment (emphasis mine)

dsolve has three extension methods, rkf45_dae, which is an extension of the rkf45 method,

I understand (and have explained) why the rkf45 method does not work, Now I do not know what the extension of the rkf45 method to rkf45_dae actually is :-(

I would have to do a lot of reading to find out. (On the other hand you could do the same reading.)

But at least you should appreciate that there is a significant difference between rkf45 and rkf45_dae and if you really want to code these explicitly then ignore rkf45 because it does not, can not and never will work for your problem. You (and possibly I) need to find out exactly what the "extensions" built-in to the rkf45_dae are.

I fixed this problem about 8hrs ago and sent a worksheet, repeated here

tp4.mw

Just in case you didn't like my original for some reason, attached to this post is your original _new_probs.mw with the same fixes in place

new_probs_fixed.mw

Care to explain what is wrong with either of these

@patient 

Sorry

I should have read your original equation more carefully. In order to apply the rkf45 method it has to be possible to convert the originall system of equations into the form

S'(x)= f₁( S(x), U(x), G(x), x)

U'(x)= f₂(S(x), U(x), G(x), x)

G'(x) =f₃(S(x), U(x), G(x), x)

with no derivatives on the right-hand side. For your system it is realtively easy to eliminate the second-order terms, but teh resulting set of equations cannot be cast in the above form. One always ends up with a derivative on the rhs. I therefore do not believe that this system of equations can be converted to a form suitable for the rkf45 method (although other methods are available)

@patient 

See the attached - is this what you want?

consys.mw

@patient 

DETools[convertsys] needs initial conditions which you don't seem to have

Typo in your fiirst equation??

sys_s:=(1/2)*x*t(diff(S(x), x))+sqrt(S(x))*U(x)*t^2+S(x)*t = 0

should be

sys_s:=(1/2)*x*t*(diff(S(x), x))+sqrt(S(x))*U(x)*t^2+S(x)*t = 0

Try using DETools[convertsys] to convert your system to first-order equations. I would have tried this but it would have meant retyping your equations (not going to happen, too time-consuming and error-prone). In future please insert them as plain text or upload a Maple worksheet using the big green arrow at the right-hand end of the toolbar icons

@momoklala 

With no changes I get

  0.00    -1.68614066179524
  0.20    -1.90456939581995
  0.40    -2.15164532806677
  0.60    -2.43167852101784

and fails betond this point

Since you were able to get solutions with S=-0.5, I tried using the contuation option to ramp this boundary condition from S=-0.5 (where it works) to S=0.5 where it currently doesn't. This produced

  0.00    -1.68614066179523
  0.20    -1.90456939581995
  0.40    -2.15164534927826
  0.60    -2.43167851062166
  0.80    -2.75018329608440
  1.00    -3.11426403117906

and failed at 1.2: which looked to me like an end-point singularity, so I then changed the method to midrich, which produced

  0.00    -1.68614066179523
  0.20    -1.90456941902650
  0.40    -2.15164539324505
  0.60    -2.43167857162100
  0.80    -2.75018331915614
  1.00    -3.11426413213177
  1.20    -3.53316438964513

For implementation, see below

tp4.mw

If you want to stick with your original worksheet, just change

S := 0.5 to S := -0.5+c and

continuation=c and method=bvp[midrich] in the dsolve options.

@spark1631 

Consensus opinion is that you have discovered a significant bug, so I have submitted a software change request. Not much help for you because history suggest that such SCRs don't tend to be dealt with very quickly :-(

@Preben Alsholm 

Must agree with Preben - again

With S=-0.5 in the worksheet I posted previously, Maple 2015 profuces

0.00    -1.18614066271840
0.20    -1.17788389906508
0.40    -1.16926705672408
0.60    -1.16025987616959
0.80    -1.15082788469627
1.20    -1.13052333151069

and Maple 18 produces

0.00    -1.18614066271840
0.20    -1.17788389906508
0.40    -1.16926705672408
0.60    -1.16025987616959
0.80    -1.15082788469627
1.20    -1.13052333151069

ie exactly the same to 13 decimal places with no sign of non-convergence. Perhaps instead of writng your own obviously non-functional code, you should use the code provided by Preben (or me), which does actually work

For this particular problem I must admit that I am slightly fascinated by the fact that the numerical answer in which you are interested is diff(diff(f(eta), eta), eta). It isn't often that one computes the numerical solution to a set of differential equations in order to obtain the numerical value of the second differential of one of the functions involved.

Still I suppose you know what you are doing - and we can only answer the questions you ask :-)

@Preben Alsholm 
I get the same answer as Preben, a slightly different way

R:= Matrix(6, 2, (i,j) -> `if`(j=1,L[i],0) ):
  AP:=NULL;
  for k to 6 do
      sol := dsolve( eval({Eq1, Eq2, bcs1, bcs2}, B = L[k]),
                     [f(eta), theta(eta)],
                     numeric,
                     output = Array([seq(j,j=0..10)]),
                     maxmesh=512,
                     AP
                   ):
      R[k,2]:=sol[2][1][1][4];
      AP:=approxsoln=sol;
  end do:
  R;


gives

0.00    -1.41421356296431
  0.20    -1.43879570418522
  0.40    -1.46322185423285
  0.60    -1.48747812977898
  0.80    -1.51155334213274
  1.20    -1.55912758973327

@Carl Love 

This problem cropped up because someone asked a question here

http://mapleprimes.com/questions/204272-Get-Rid-Of-Absolute-Value-When-Solve

looking to use solve() rather naively for what seemed to be a fairly harmless-looking set of equations.

sys := {x = abs(y-1), 5^(2*x+2) = 25^(2*x-3), 5^(3*x+1)/5^(x-1) = 25^(2*x-3)};

As it happens this system is trivial to solve: the second and third equations are in fact identical and contain only the variable x, so one can get x (=4)  by either solve() or isolate() on sys[2], then substitute in sys[1] for y (=5,-3).

But for something this simple, I was surprised/disappointed that solve(sys), returned

Error, (in testeq) numeric exception: overflow

It was whilst digging around into why this might happen that I came up with the equation I originally posted - one can get it by

simplify(subs(sys[1],sys[2]));

OP was not too impressed at the performance of Maple's solve() command on this system - and I had to agree with him

@smith_alpha 

You have to ask yourself a serious question about working methods.

1. One choice is to create a utility file which can be read/executed by multiple worksheets and edited by none of them
2. The other choice is to create a utility file which can be read/executed by multiple worksheets, and edited by one of them (aka the magic worsheet)

Either works: I prefer the first. Probably because I would never remember which was the *magic* worksheet with edit privilege. If I did come up with an improvement to the utility whilst working on ws99.mw, would I edit it as a text file applied across all worksheets, or remember which of the hundred or so worksheets happened to define it snd edit it there??

Essentially it is a question of how you handle the primary reference for a utility: is it the utility itself, or a worksheet which defines that utility. Feel free to pick one