# Getting error message "The minimal damping factor of 1/10000 has been reached"

As a very beginner into Mathematica, I encounter error message as shown below when I am solving the system of three second order ODEs numerically using ParametricNDSolveValue:

FindRoot::dfmin: The minimal damping factor of 1/10000 has been reached.

ParametricNDSolveValue::fempsf: PDESolve could not find a solution.

Below is the essential part of my codes, starting with list of constants,

(*List of constants*)
w = 1;
p = 1;
q = 0.05;
u = 1;
Rba = 0;
Rca = 1;
Kb = 1.38*10^-23;
T = 25 + 273;
miumaxC = 0;
alphaA = 20;
alphaB = 20;
alphaC = 1;
betaA = 1;
betaB = 1;
betaC = 1;
kgA = 1;
kgB = 1;
kgC = 1;


Intermediate functions,

(*intermediate functions*)
cp1[x_] = miumaxA/2*(1 - Tanh[alphaA (x - betaA)])*Kb*T;
cp2[x_] = miumaxB/2*(1 - Tanh[alphaB (x - betaB)])*Kb*T;
cp3[x_] = miumaxC/2*(1 - Tanh[alphaC (x - betaC)])*Kb*T;
diffintermediate1[x_] = a[x]*x^2*cp1'[x];
diffintermediate2[x_] = b[x]*x^2*cp2'[x];
diffintermediate3[x_] = (phi[x] - 1)*x^2*cp3'[x];


Defining ODEs and boundary conditions before solving,

(*Dimensionless species balance equations with neumann BCs for \
species A,B,and C*)
ODE1 = a''[x] + 2/x*a'[x] + 1/(Kb*T)/x^2*diffintermediate1'[x] -
TM^2*phi[x]*a[x] ==
NeumannValue[0, x == 0] +
NeumannValue[kgA*(a[x] - 1), x == 1]; ODE2  =
b''[x] + 2/x*b'[x] + 1/(Kb*T)/x^2*diffintermediate2'[x] +
TM^2*phi[x]/p*(a[x] - w*b[x]) ==
NeumannValue[0, x == 0] + NeumannValue[kgB/p* (b[x]), x == 1];
ODE3  = q*(phi''[x] + 2/x*phi'[x] +
1/(Kb*T)/x^2*diffintermediate3'[x]) - TM^2*u*w*phi[x]*b[x] ==
NeumannValue[0, x == 0] + NeumannValue[kgC/q* (phi[x] - 1), x == 1];

(*Grouping the ODEs for input in NDSolveValue function*)
ODEs = {ODE1, ODE2 , ODE3};

(*Solving the system of ODES*)
solutions =
ParametricNDSolveValue[{ODEs}, {a, b, phi}, {x, 0, 1}, {TM,
miumaxA, miumaxB}];


I would get the abovementioned error message when I try to calculate and plot the solutions at TMplot = 10 to investigate the effect of different miumax. Below is the code for plotting the solutions:

MiumaxList = {-4,-2,0,2,4};
(*Solution profile plotting at different miumaxA*)
TMplot  =10;

(*Creating legend colours*)
legendColours = {};
legendEntries = {};
For[i=1,i<=Length[MiumaxList],i++,
legendColours = Insert[legendColours,RGBColor[0.2*i,0.35,0.7],i];
]

(*Creating list of solutions for varying Miumax*)
solutionListA = {};
solutionListB = {};
solutionListPhi = {};
For[i=1,i<=Length[MiumaxList],i++,
solutionListA= Insert[solutionListA,solutions[TMplot,MiumaxList[[i]],MiumaxList[[i]]][[1]][x] ,i];
solutionListB = Insert[solutionListB,solutions[TMplot,MiumaxList[[i]],MiumaxList[[i]]][[2]][x] ,i];
solutionListPhi = Insert[solutionListPhi,solutions[TMplot,MiumaxList[[i]],MiumaxList[[i]]][[3]][x] ,i];
];

(*Plotting of solutions a, b, and phi for varying Miumax*)
Plot[solutionListA,{x,0,1},PlotRange->{Automatic,Full},PlotLabel->"a (x) at different miumax pair, TM = "[TMplot],AxesLabel->{"x","a(x)"},PlotStyle->legendColours,PlotLegends->LineLegend[legendColours,MiumaxList]]
Plot[solutionListB,{x,0,1},PlotRange->{Automatic,Full},PlotLabel->"b (x) at different miumax pair, TM = "[TMplot],AxesLabel->{"x","a(x)"},PlotStyle->legendColours,PlotLegends->LineLegend[legendColours,MiumaxList]]
Plot[solutionListPhi,{x,0,1},PlotRange->{Automatic,Full},PlotLabel->"a (x) at different miumax pair, TM = "[TMplot],AxesLabel->{"x","a(x)"},PlotStyle->legendColours,PlotLegends->LineLegend[legendColours,MiumaxList]]


There are two things you can try. In version 13.1 parametric functions when used with the finite element method can take a new initial seeding. This is relevant in your case as the nonliterary is very strong as numbers increase and a good initial seed might help to push the solutions further out.

Let's have a look:

(*Solving the system of ODES*)
solutions =
ParametricNDSolveValue[{ODEs}, {a, b, phi}, {x, 0, 1}, {TM, miumaxA,
miumaxB}];


We can compute a solution for this set of parameters:

oldSolution = solutions[3, -4, -4]
(* InterpolatingFunction[.... *)


But when we increase the 3 to a 4 we get a convergence error:

solutions[4, -4, -4]


What you can do now is to use the oldSolution as an initial seed for the next solution like so:

solutions[4, -4, -4,
"InitialSeeding" ->


I can not say how far you can push these out but give it a shot. If this does not work, not all hope is lost. There is a second way to do this. In that case you convert your stationary ODEs to time dependent ODEs and time integrate them for a pseudo time from 0 to 1 and use for example the first parameter like

NDSolve[ tsolutions[t*10, -4, -4],..., {t,0,1},...]


This has the advantage that the adaptive time stepping algorithm of NDSolve will take care of the step size selection; on the downside you need to set up a time dependent problem and initial conditions.

Both approaches are shown in the monograph on Hyperelasticity.

Good luck!

• Thank you so much for the advice. I just tried the first method, which is to alter the initial seeding by calling it using solutions[4, -4, -4, "InitialSeeding" -> Thread[Through[{a, b, phi}[x]] == Through[oldSolution[x]]]]. However, I would get error: "ParametricNDSolveValue::fpct: Too many parameters in {TM,miumaxA,miumaxB} to be filled from ... ..." It seems like I have to change the initial seedings in ParametricNDSolveValue[... ...]. Is there any workaround that changes the initial seeding after declaring ParametricNDSolveValue? Commented Jul 11, 2022 at 9:48
• @joe, what is you \$Version? Commented Jul 11, 2022 at 10:48
• I am currently using ver 13.0.1.0 Commented Jul 11, 2022 at 12:22
• @joe, as mentioned in the beginning of the post, the above works in version 13.1. In version 13.0 you'd have to try the second method that I outlined. Make this a time dependent problem and add a factor of t in front of the parameters. Commented Jul 11, 2022 at 13:13