# Solving stiff equations in Mathematica

I have problem to solve stiff equations. Any idea on how to solve this? I have tried "StiffSwitching" but it didnt work. Im solving this using Mathematica 10.

Here is my code. Im sorry if I wrote the code badly but I have fllowed the instructions on how to copy and paste code from Mathematica and this is what I end up with.

With[{n = 0.6}, ODE3[L_] := {Derivative[1][H][\[Eta]] == -2*F[\[Eta]] - ((1 - n)/(1 + n))*\[Eta]*Derivative[1][F][\[Eta]],
F[\[Eta]]^2 - (G[\[Eta]] + 1)^2 + (H[\[Eta]] + ((1 - n)/(1 + n))*\[Eta]*F[\[Eta]])*Derivative[1][F][\[Eta]] ==
(Derivative[1][F][\[Eta]]^2 + Derivative[1][G][\[Eta]]^2)^((n - 1)/2)*Derivative[2][F][\[Eta]] +
Derivative[1][F][\[Eta]]*((n - 1)*(Derivative[1][F][\[Eta]]^2 + Derivative[1][G][\[Eta]]^2)^((n - 3)/2)*(Derivative[1][F][\[Eta]]*Derivative[2][F][\[Eta]] +
Derivative[1][G][\[Eta]]*Derivative[2][G][\[Eta]])), 2*F[\[Eta]]*(G[\[Eta]] + 1) + (H[\[Eta]] + ((1 - n)/(1 + n))*\[Eta]*F[\[Eta]])*Derivative[1][G][\[Eta]] ==
(Derivative[1][F][\[Eta]]^2 + Derivative[1][G][\[Eta]]^2)^((n - 1)/2)*Derivative[2][G][\[Eta]] +
Derivative[1][G][\[Eta]]*((n - 1)*(Derivative[1][F][\[Eta]]^2 + Derivative[1][G][\[Eta]]^2)^((n - 3)/2)*(Derivative[1][F][\[Eta]]*Derivative[2][F][\[Eta]] +
Derivative[1][G][\[Eta]]*Derivative[2][G][\[Eta]])), F[0] == 0, G[0] == 0, H[0] == 0, F[L] == 0, G[L] == 0}]

sol4[L_] := Flatten[NDSolve[ODE3[L], {F, G, H}, \[Eta],
Method -> {"Shooting", "StartingInitialConditions" -> {G[0] == 0, Derivative[1][G][0] == -0.65, F[0] == 0, Derivative[1][F][0] == 0.55,
H[0] == 0}}]]

Plot[Evaluate[{F[\[Eta]], G[\[Eta]], H[\[Eta]]} /. sol4[14]], {\[Eta],
0, 14}, PlotStyle -> Thick, Frame -> True,
FrameLabel -> {"\[Eta]", "F, G, H"},
Epilog -> {Dashed, Thickness[0.005],
Line[{{0, 0.8844}, {8, 0.8844}}]}]


I got this type of error:

At \[Eta] == 11.976704516552642, step size is effectively zero; singularity or stiff system suspected.


This only works if I let n=1. But I have to plot for different values of n ranging from 0.6 to 1.4. Below is what plots I got from Mathematica.

It would be great if someone could help me and explain to me. Thank you in advance!

• Have you checked if the solution of this BVP is singular there? Can you at least show the plots? It's not completely clear that it is really stiffness that causes this. If telling Mathematica to switch to a stiff solver didn't help, the problem could be elsewhere. – Kirill Jan 27 '16 at 23:16
• Hi. I already edited my question above. I have added more of the coding and also the plot which I got. About the singularity, I am not really sure about that. How do I know if this BVP is singular btw? Thank you. – dayana nisa Jan 27 '16 at 23:33

When solving a boundary value problem, the shooting method uses nonlinear rootfinding to find the initial conditions that would satisfy the boundary conditions. To see what's going on here, replace your BVP with the IVP corresponding to the first guess that you told the shooting method to start with. That is, replace

F[L] == 0, G[L] == 0


with

F'[0] == 0.55, G'[0] == -0.65


and plot the resulting solutions of the IVP. Unlike when solving the BVP, Mathematica will actually give you the solution up to time 11.977 when the singularity occurs. You can clearly see in the plot how the functions blow up — so it is definitely not stiffness causing this:

This makes it clear that your problem is that certain (many?) initial values result in singularities. So the shooting method is not suitable, unless by some coincidence the initial values that it tries just happen to not lead to singularities.

When I tried setting the starting initial conditions to F'[0] == 0.5, G'[0] == -0.7` (a bit of guesswork and trial and error), it seemed to work correctly, giving a solution that satisfies the BCs: