# How to solve a ParametricNDSolve issue when the system seems to be stiff?

I have a problem with ParametricNDSolve[]. Mathematica gives the warning that a stiff system is suspected. I have varied the methods but unfortunately I was not able to solve it.

First some formulas that are used in the definition of the differential equation within the ParametricNDSolve[] statement:

r = 0.005; η = 1;
P[t_, T0_] := Exp[-r (T0 - t)]
Bx[x_, α_, γ_] := α/γ (1 - Exp[-γ x])
T1 = T0 + 1/365;


This is a list of test parameters:

testPars = {α -> 0.2563, γ ->
0.3753, ρ13 -> -0.2414, ρ12 -> -0.7369, ρ23 ->
0.1300, σv -> 2.5618, σS -> 0.6143, κ ->
5.7734, v0 -> 0.1559};


and finally the ParametricNDSolve[] statement as how I have entered it:

nfun = n /. ParametricNDSolve[{\!$$\*SubscriptBox[\(∂$$, $$τ$$]\ $$n[τ, u]$$\) == n[τ, u] (-κ + u σv (ρ13 σS + ρ23 Bx[T1 - T0 + τ, α, γ])) + 1/2 n[τ, u]^2 σv^2 + 1/2 (u^2 - u) (σS^2 + Bx[T1 - T0 + τ, α, γ]^2 + 2 ρ12 σS Bx[ T1 - T0 + τ, α, γ]), n[0, u] == 0}, n, {τ, 0, 5}, {u, 0, 10^4}, {α, γ, ρ13, ρ12, ρ23, σv, σS, κ}, Method -> "StiffnessSwitching"]


I get an error when I plot the solution for the testparameter set:

Plot[nfun[α, γ, ρ13, ρ12, ρ23, σv, σS, κ][τ, 0] /. testPars, {τ, 0, 1}]


Does anyone have an idea on how to solve it? Your help is much appreciated.

-
Thanks for letting me know Nasser, I have adjusted my post. –  Frits Feb 3 at 14:36
You may consider exploring your ODE with a smaller domain for u. For example when I use ParametricNDSolve[…,{t,0,5},{u,0,5}], I don't see any issue with the stiffness. Your solution may also be better visualized with Plot3D. –  leibs Feb 4 at 5:49
After a little more poking around it looks like as u approaches 7 the ODE blows up to infinity. As this isn't an issue with Mathematica you may consider bringing this up over at math.stackexchange.com for further discussion. They might be able to help classify the ODE and figure out if this is the behavior you should be seeing. –  leibs Feb 4 at 5:57