# How to solve stiff forth order ODE using NDSolve

I'm trying to get solution of the following ode using following code. But no plotting appears. I need help to solve the problem.

a = 0.5;
b = 2.0;
x = 0.1;
t = 0.1;
Ha = 1.0;
Q0 = 1.0;
we = 0.2;
alp = 1.0;
F = Q0*Exp[-alp*t]
dVals = {0.01, 0.1, 0.2, 0.3};
h = 1.0 - (a*Cos[Pi*(x - t)]*Cos[Pi*(x - t)]);
sys2 = {(1.0 - d) s''''[y] + (2.0*d*we*s''[y]*s''''[y]) + (2.0*d*we*
s''[y]*s''[y]) == Ha*Ha*s'[y]};
iv2 = {s[-h] == -F/2.0, s[h] == F/2.0, s'[-h] == 0.0, s'[h] == 0.0};
sol2[d_] :=
NDSolve[Join[sys2, iv2], s, {y, -h, h},
Method -> {"StiffnessSwitching",
Method -> {"ExplicitRungekutta", Automatic}},
WorkingPrecision -> 35, AccuracyGoal -> 5, PrecisionGoal -> 4];
Plot[Evaluate[Table[s'[y] /. sol2[d], {d, dVals}]], {y, -1, 1},
PlotRange -> Full, AxesLabel -> {y, \[CapitalPsi]}, Frame -> False,
AxesStyle -> Black,
PlotStyle -> {{Line, Black}, {DotDashed, Black}, {Dotted,
Black}, {Dashed, Black}},
PlotLegends ->
Placed[{"d = 0.01  ", "d = 0.1  ", "d = 0.2  ", "d = 0.3 "},
Center]]

• we is undefined. Commented Jul 19, 2017 at 12:47
• Value of we is very small. It may be take as 0.2. I edited the question. Hope it will be okay to get a solution now. Commented Jul 19, 2017 at 13:19

It is good practice to use rational rather than decimal numbers with NDSolve when possible. Also, the equation actually is not stiff, despite the error message that sometimes occurs. So, use

a = 1/2; b = 2; x = 1/10; t = 1/10; Ha = 1; Q0 = 1; alp = 1; F = Q0*Exp[-alp*t]; we = 2/10;
dVals = {1/100, 1/10, 2/10, 3/10};
h = 1 - (a*Cos[Pi*(x - t)]*Cos[Pi*(x - t)]);
sol2[d_] := Flatten@NDSolve[{(1 - d) s''''[y] + (2*d*we*s''[y]*s''''[y]) + (2*d*
we*s''[y]*s''[y]) == Ha*Ha*s'[y], s[-h] == -F/2, s[h] == F/2, s'[-h] == 0, s'[h] == 0},
{s, s', s'', s'''}, {y, -h, h}, WorkingPrecision -> 30]
Plot[Evaluate[Table[s'[y] /. sol2[d], {d, dVals}]], {y, -h, h},
PlotRange -> Full, AxesLabel -> {y, Ψ}, Frame -> False, AxesStyle -> Black,
PlotStyle -> {{Line, Black}, {DotDashed, Black}, {Dotted, Black}, {Dashed, Black}},
PlotLegends -> Placed[{"d = 0.01 ", "d = 0.1 ", "d = 0.2 ", "d = 0.3 "}, Center]]


Note, also, that the Plot range has been reduced to correspond to the NDSolve range.

It may seem strange that s'[t] is essentially unchanged as d is increased from near zero to 3/10. In fact, s[t] and s''[t] also are essentially unchanged. s'''[t], on the other hand changes significantly, as shown here.

Plot[Evaluate[Table[s'''[y] /. sol2[d], {d, dVals}]], {y, -h, h},
PlotRange -> Full, AxesLabel -> {y, Ψ}, Frame -> False, AxesStyle -> Black,
PlotStyle -> {{Line, Black}, {DotDashed, Black}, {Dotted, Black}, {Dashed, Black}},
PlotLegends -> Placed[{"d = 0.01 ", "d = 0.1 ", "d = 0.2 ", "d = 0.3 "}, Center]]


The third derivative is singular for y near h for d == 3/10, because the coefficient of s''''[y] vanishes there,

Coefficient[(1 - d) s''''[y] + (2*d*we*s''[y]*s''''[y]) + (2*d*we*
s''[y]*s''[y]) - Ha*Ha*s'[y], s''''[y]]
(* 1 - d + 2/5 d s''[y] *)


Indeed, the coefficient vanishes for d == 3/10 and s''[h] approximately equal to -5.4. Empirically, the solution is on the verge of becoming singular at y == h for d == 2985209/10000000.

Plot[Evaluate[{s[y], s'[y], s''[y], s'''[y] + 10} /.   sol2[2985209/10000000]],
{y, -h, h}, PlotLegends -> Placed[{s[y], s'[y], s''[y], s'''[y] + 10}, {.8, .8}]]


• Thank you. I have to show how the solution is affected by the two key parameters 'd' and 'we' in two different plots. In this case no charges take place in the above plot. Then what is the advantages of taking such parameter in ode. I'm really confused. Commented Jul 19, 2017 at 13:38
• @BiswajitMallick With larger values of we, the answer is likely to change, although solving the ODE will be more difficult. Commented Jul 19, 2017 at 13:54
• @BiswajitMallick I'll look more at this. Commented Jul 19, 2017 at 15:05
• Thank you very much. @bbgodfrey Commented Jul 19, 2017 at 16:03