# Unable to set MaxStepSize when boundary conditions are specified using NeumannValue and DirichletCondition

Short story: I am trying to specify the maximum step size for NDSolve function, but it appears to be working only if I do not specify boundary conditions using NeumannValue and DirichletCondition functions.

Long story: I have the following set of equations and boundary conditions

Gc0 = DirichletCondition[c[y] == cB, y == -1];
GmA = NeumannValue[- (Q[y] + (2 + c[y])/3), y == -1];
GmB = NeumannValue[(Q[y] + (2 + c[y])/3), y == 1];
GQA = NeumannValue[-(4/15) m[y], y == -1];
GQB = NeumannValue[(4/15) m[y], y == 1];

eqs1 = {
- c'[y] + m[y] == 0, Gc0,
- m''[y] + Q'[y] + c'[y]/3 + m[y] == GmA + GmB,
- Q''[y] + (4/15) m'[y] + Q[y] == GQA + GQB
}


So I set up the solver i the standard way, with the max step size being 2/100

nSol[eq_, par_] := NDSolve[eq /. par, {c, m, Q}, {y, -1, 1}, MaxStepSize -> 2/100]
res = nSol[eqs1, {cB -> 0.03}];
Plot[{c[y] /. res, m[y] /. res, Q[y] /. res}, {y, -1, 1}, PlotRange -> All]


The above works completely fine, but when I checked the grid of the interpolating function I find

c["Grid"] /. res

{{{-1.}, {-0.9}, {-0.8}, {-0.7}, {-0.6}, {-0.5}, {-0.4}, {-0.3},{-0.2}, {-0.1}, {0.}, {0.1}, {0.2}, {0.3}, {0.4}, {0.5}, {0.6},{0.7}, {0.8}, {0.9}, {1.}, {-0.95}, {-0.85}, {-0.75}, {-0.65}, {-0.55}, {-0.45}, {-0.35}, {-0.25}, {-0.15}, {-0.05}, {0.05}, {0.15}, {0.25}, {0.35}, {0.45}, {0.55}, {0.65}, {0.75}, {0.85}, {0.95}}}


Which does not look like spacing from -1 to 1 with the max step size being 2/100.

However, if I set my equations in the following way

eqs2 = {
- c'[y] + m[y] == 0,
- m''[y] + Q'[y] + c'[y]/3 + m[y] == 0,
- Q''[y] + (4/15)m'[y] + Q[y] == 0,

c[-1] - cB == 0,
m'[-1] - Q[-1] -(2+c[-1])/3 == 0,
m'[1] - Q[1] -(2+c[1])/3 == 0,
Q'[-1] - (4/15)m[-1] == 0,
Q'[1] - (4/15)m[1] == 0
}


and I repeat the calculation I get the same plots, but the MaxSetpSize now seem to be doing what I would expect

c["Grid"] /. res

{{{-1.}, {-0.999904}, {-0.999808}, {-0.995684}, {-0.99156}, {-0.987436}, {-0.974794}, {-0.962152}, {-0.94951}, {-0.936868}, {-0.916868}, {-0.896868}, {-0.876868}, {-0.856868}, {-0.836868}, {-0.816868}, {-0.796868}, {-0.776868}, {-0.756868}, {-0.736868}, {-0.716868}, {-0.696868}, {-0.676868}, {-0.656868}, {-0.636868}, {-0.616868}, {-0.596868}, {-0.576868}, {-0.556868}, {-0.536868}, {-0.516868}, {-0.496868}, {-0.476868}, {-0.456868}, {-0.436868}, {-0.416868}, {-0.396868}, {-0.376868}, {-0.356868}, {-0.336868}, {-0.316868}, {-0.296868}, {-0.276868}, {-0.256868}, {-0.236868}, {-0.216868}, {-0.196868}, {-0.176868}, {-0.156868}, {-0.136868}, {-0.116868}, {-0.0968684}, {-0.0768684}, {-0.0568684}, {-0.0368684}, {-0.0168684}, {0.00313162}, {0.0231316}, {0.0431316}, {0.0631316}, {0.0831316}, {0.103132}, {0.123132}, {0.143132}, {0.163132}, {0.183132}, {0.203132}, {0.223132}, {0.243132}, {0.263132}, {0.283132}, {0.303132}, {0.323132}, {0.343132}, {0.363132}, {0.383132}, {0.403132}, {0.423132}, {0.443132}, {0.463132}, {0.483132}, {0.503132}, {0.523132}, {0.543132}, {0.563132}, {0.583132}, {0.603132}, {0.623132}, {0.643132}, {0.663132}, {0.683132}, {0.703132}, {0.723132}, {0.743132}, {0.763132}, {0.783132}, {0.803132}, {0.823132}, {0.843132}, {0.863132}, {0.883132}, {0.903132}, {0.923132}, {0.943132}, {0.963132}, {0.981566}, {1.}}}


I also tried different values of MaxStepSize, but it seems to be doing something only if I define the boundary conditions as in eqs2 and not as in eqs1.

What I am missing here?

• The DirichletCondition and NeumannValue trigger the use of the finite element method, which leads to a spatial discretization of your ode. Now, MaxStepSize is an option for a temporal integrator in NDSolve. If you want to time integrate the ODE (this is probably what you want) then you should rewrite the DirichletCondition and NeumannValue in equational form. Aug 11, 2021 at 12:22