# Error test failure when solving two coupled ODEs

by NDSolve; this is work until y=6.22

sol = NDSolve[{f'''[y] + f[y] f''[y] + 4 - (f'[y])^2 == 0,
g''[y] + 0.01*f[y] g'[y] == 0, f[0] == 0, f'[0] == 0, g'[0] == -1,
f'[10] == 2, g[10] == 0}, {f, g}, {y, 0, 10}]
g[0] /. sol


after that I got this

Error test failure at y == 6.0045645472974165; unable to continue.

I need to increase the value ...is there another way rather than ND-solve.

• i cannot increase the max value here ..bvb's involving ode's – One Two Dec 11 '19 at 19:15
• People here generally like users to post code as Mathematica code instead of just images or TeX or regular math notation, so they can copy-paste it. It makes it convenient for them and more likely you will get someone to help you. It also helps to post working code (that is, the complete NDSolve command you used in this instance). You may find this meta Q&A helpful. – Michael E2 Dec 11 '19 at 19:25
• thank you so much – One Two Dec 11 '19 at 20:00
• Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. – bbgodfrey Dec 12 '19 at 14:01

This ODE system can be solved using the option Method -> "StiffnessSwitching", although the computation proceeds slowly.

sol = NDSolveValue[{f'''[y] + f[y] f''[y] + 4 - (f'[y])^2 == 0,
g''[y] + 0.01*f[y] g'[y] == 0, f[0] == 0, f'[0] == 0, g'[0] == -1,
f'[10] == 2, g[10] == 0}, {f[y], g[y]}, {y, 0, 10},
Method -> "StiffnessSwitching", MaxSteps -> 10^6];
Plot[sol, {y, 0, 10}, ImageSize -> Large, AxesLabel -> {y, "f,g"},
LabelStyle -> {Black, Bold, 15}]


Particular characteristics of the ODE system permit solving it more simply than with the general approach given above. Specifically, g[y] does not appear in either of the ODEs, although its derivatives do. Consequently, the boundary condition, g[10] == 0 can be satisfied merely by solving the system with g[0] == 0 satisfied instead of g[10] == 0, and then subtracting from the solution the value of g[10] as just calculated to force g[10] == 0. To proceed, solve the system parametrically with

s = ParametricNDSolveValue[{f'''[y] + f[y]*f''[y] + 4 - (f'[y])^2 == 0, f[0] == 0,
f'[0] == 0, f''[0] == fpp, g''[y] + coef*f[y]* g'[y] == 0, g[0] == 0, g'[0] == -1},
{f[y], g[y] - g[10], f'[10]}, {y, 0, 10}, {fpp, coef},
Method -> "StiffnessSwitching", WorkingPrecision -> 30];


A coupling coefficient, coef, has been introduced to simplify obtaining solutions for various coupling coefficients, as requested by the OP in a comment below. (WorkingPrecision -> 30 is needed only for values of coef of order 4 or larger.) Note also that ParametricNDSolve is instructed to return g[y] - g[10] as explained just above, and also f'[10]. The desired solution is obtained by varying fpp with FindRoot until f'[10] == 2.

First, plot f'[10] as a function of fpp.

Plot[Last[s[fpp0, 1/100]] - 2, {fpp0, 3.36, 3.5}, ImageSize -> Large,
AxesLabel -> {fpp, "f'[10]"}, LabelStyle -> {Black, Bold, 15}]


Perhaps surprisingly, there are two solutions for f'[10] == 2, at fpp of about 3.37 and 3.49. The first of these yields the solution

fpp0 /. FindRoot[Last[s[fpp0, 1/100]] == 2., {fpp0, 3.37}, Evaluated -> False]
(* 3.37136 *)
Plot[Evaluate@Most@s[%, 1/100], {y, 0, 10}, ImageSize -> Large,
AxesLabel -> {y, "f,g"}, LabelStyle -> {Black, Bold, 15}]


producing a plot identical to that above, as expected. The new solution is

fpp0 /. FindRoot[Last[s[fpp0, 1/100]] == 2., {fpp0, 3.49}, Evaluated -> False]
(* 3.48628 *)
Plot[Evaluate@Most@s[%, 1/100], {y, 0, 10}, ImageSize -> Large,
AxesLabel -> {y, "f,g"}, LabelStyle -> {Black, Bold, 15}]


I have obtain solutions for coef as large as 7, shown below, without difficulty. Still larger values require a higher WorkingPrecision.

• thank you so much ...your help will change many things related to my success – One Two Dec 12 '19 at 10:09
• but still cannot give me the value of g[0]/.sol – One Two Dec 12 '19 at 10:31
• Happy to help. Use Last@sol /. y -> 0, which yields 8.95717. – bbgodfrey Dec 12 '19 at 14:01
• in this code i faced a lot pf problem when I replace 0.01 by number greater than 5 or numbers like 0.72 ... can i solve this system by another built in function ....sorry for keep asking .. – One Two Dec 13 '19 at 21:52
• @OneTwo A modification of my earlier answer will handle the larger coupling constant. I shall add it in a few hours. – bbgodfrey Dec 13 '19 at 23:33

We can simplify the procedure since the f-equation can be solved totaly independent of the g-function. Thanks to the inspiration by @bbgodfrey.

First solve the f-equation with initial conditions. (I found parameters i1 and i2 with FindRoot since shooting method did not work with my version 8.0 in that case.)

fsol[i_?NumericQ] :=
f /. First@
NDSolve[{f'''[y] + f[y] f''[y] + 4 - (f'[y])^2 == 0, f[0] == 0,
f'[0] == 0, f''[0] == i}, f, {y, 0, 10}]

i1 = i /. FindRoot[Derivative[1][fsol[i]][10] == 2, {i, 3}]

(*   3.37136   *)

i2 = i /. FindRoot[Derivative[1][fsol[i]][10] == 2, {i, 5}]

(*   3.48628   *)

Plot[Evaluate[{fsol[i1][y], fsol[i2][y]}], {y, 0, 10},
PlotStyle -> {Green, Red}, PlotRange -> All]


Now use the found fsol[y] interpolating function for the second g-equation together with the cuppling constant coef all for the both found i1 and i2.

gsol[i_?NumericQ, coef_?NumericQ] :=
g /. First@
NDSolve[{g''[y] + coef*fsol[i][y]* g'[y] == 0, g'[0] == -1,
g[10] == 0}, g, {y, 0, 10}]

gsol[i1, .01]

Manipulate[
Plot[Evaluate[{gsol[i1, coef][y], gsol[i2, coef][y]}], {y, 0, 10},
PlotStyle -> {Green, Red}, PlotRange -> All], {{coef, .01}, 0, 10}]


Edit

Appendix to the OP's comment ..that to see the influence of changing the parameter( Coeff of f(y)g'(y)) from 0.01 to 0.02 .......4 5 10 100 1000 and so on . so i need numerical result for g(0) f''(0) with different value of the parameter .. .

(Needs higher workingPrecision for fsol and rationalizing:)

fsol[i_?NumericQ] :=
f /. First@
NDSolve[{f'''[y] + f[y] f''[y] + 4 - (f'[y])^2 == 0, f[0] == 0,
f'[0] == 0, f''[0] == Rationalize[i, 0]}, f, {y, 0, 10},
WorkingPrecision -> 25]

(vals = Table[{coef,
Evaluate[{gsol[i1, coef][0]*Derivative[2][fsol[i1]][0],
gsol[i2, coef][0]*Derivative[2][fsol[i2]][0]}]}, {coef, {0.01,
0.02, .05, .1, .2, .5, 1, 2, 5, 10, 100, 1000}}]) // TableForm


Plot[Evaluate[{gsol[i1, coef][0]*Derivative[2][fsol[i1]][0],
gsol[i2, coef][0]*Derivative[2][fsol[i2]][0]}], {coef, 1/100,
1000}, PlotStyle -> {Green, Red}]
`
• thank you so much for your help.... the aim from solving this system of equations ..that to see the influence of changing the parameter( Coeff of f(y)g'(y)) from 0.01 to 0.02 .......4 5 10 100 1000 and so on . so i need numerical result for g(0) f''(0) with different value of the parameter ....what ever you help a lot and will be selfish to ask more about this ....again thank you so much – One Two Dec 14 '19 at 7:21
• I made an edit to take this into account. – Akku14 Dec 14 '19 at 8:36
• again and again thank to all .... – One Two Dec 14 '19 at 10:33
• but why this value Derivative[2][fsol[i2]][0] dose not change as the parameter change ...and by the way i need the values for g(0),..f''(0)..separately – One Two Dec 14 '19 at 10:35
• Think a little bit! – Akku14 Dec 14 '19 at 10:39