# NDSolve Convergence test failure and significant effect of DifferenceOrder on eventual results

I am solving a non linear partial differential equation with what I call free boundary conditions (solid mechanicists would know this as simply supported).

I realized that this boundary condition does induce a certain stiffness in my equation which is not so strong for, say, periodic boundary condition. How do I know this? The simulation runs longer for free boundaries with a stiff method such as LSODA. This tells me that the time step is being made smaller and smaller to accommodate this rapidly changing time variation.

I also have various structures forming as obvious from the plot of my result. So I use the TensorProductGrid spatial discretization method as I gather it would be useful to solve this equation with a preferential spatial grid treatment (denser grid in regions of necessity).

As I increase the DifferenceOrder for my method/spatial discretization, I notice that I run into convergence issues NDSolve::ndcf at a later stage. Viz., for a smaller DifferenceOrder I have a convergence issue at a later time. This initially struck me as being a little strange as wouldn't a superior difference order help alleviate convergence? But then obviously, a greater difference order would catch the issue with convergence earlier!

## Here is my code:

m=0.05;
\[Epsilon]=10^-6;
Bo=1/300;
\[Delta]=10^-7;
Bi=1;
K1=1;
wn=2;

(*R=Ra/(Pr S1);*)
R=0*1/10;
{xMin,xMax}=wn{0,79.5788};
k=0.0677*92.389/(xMax/wn);
TMax=10000*100;
n=1;

uSolpbc[t_,x_]=u[t,x]/.NDSolve[{\!$$\*SubscriptBox[\(\[PartialD]$$, $$t$$]$$u[t, x]$$\)+\!$$\*SubscriptBox[\(\[PartialD]$$, $$x$$]$$( \*SuperscriptBox[\(u[t, x]$$, $$3$$]\
\*SubscriptBox[$$\[PartialD]$$, $$x, x, x$$]u[t, x])\)\)-Bo \!$$\*SubscriptBox[\(\[PartialD]$$, $$x$$]$$( \*SuperscriptBox[\(u[t, x]$$, $$3$$]\
\*SubscriptBox[$$\[PartialD]$$, $$x$$]u[t, x])\)\)+m  \!$$\*SubscriptBox[\(\[PartialD]$$, $$x$$]$$( \*SuperscriptBox[\(( \*FractionBox[\(u[t, x]$$, $$K1 + Bi\ u[t, x]$$])\), $$2$$]\
\*SubscriptBox[$$\[PartialD]$$, $$x$$]u[t, x])\)\) +\[Epsilon]/(Bi u[t,x] + K1) +\[Delta] u[t,x]^3/(Bi u[t,x] + K1)^3 D[D[u[t,x],x],x] + R \!$$\*SubscriptBox[\(\[PartialD]$$, $$x$$]$$( \*SuperscriptBox[\(u[t, x]$$, $$3$$]\
\*SubscriptBox[$$\[PartialD]$$, $$x$$]u[t, x])\)\)== 0,

u[0,x]==1-0.05 (Cos[k*x]),

(*PERIODIC BOUNDARY CONDITIONS-START*)
(*u[t,xMin]== u[t,xMax],
Derivative[0,1]u[t,xMin]==Derivative[0,1]u[t,xMax],
Derivative[0,2]u[t,xMin]==Derivative[0,2]u[t,xMax],
Derivative[0,3]u[t,xMin]==Derivative[0,3]u[t,xMax]*)
(*PERIODIC BOUNDARY CONDITIONS-END. Uncomment to run periodic case*)

(*FREE BOUNDARY CONDITIONS-START*)

Derivative[0,1][u][t,xMin]==0,Derivative[0,1][u][t,xMax]==0,Derivative[0,3][u][t,xMin]==0,Derivative[0,3][u][t,xMax]==0
(*FREE BOUNDARY CONDITIONS-START*)

},
u,
{t,0,TMax},
{x,xMin,xMax},
(*PrecisionGoal->Automatic,*)
Method->{"MethodOfLines","Method"->"LSODA","TemporalVariable"->t,"SpatialDiscretization"->{"TensorProductGrid","MinPoints"->10000,"MaxPoints"->10000,"DifferenceOrder"->1}}][[1]]


## And the plotting function:

evapfbc2 = Plot[
uSolpbc[6981*100, x], {x, xMin, xMax},
PlotStyle -> {Black, Thick, Dashed},
AxesLabel -> {"\[Lambda]", "h"},
BaseStyle -> {FontWeight -> "Plain", FontSize -> 18},
PlotRange -> {{0, xMax}, {0, 3.5}}
]


## My questions are:

1. Should I be increasing the difference order for better results? (I understand that with a difference order set to n, it allows the stiff solver "elbow room" up to and including order n)
2. Is this the wrong method to use and is there something I could do in mathematica that would make this better? I am assuming that I would need to use a stiff solver for this problem.

## Note:

With a difference order of 1, this runs for about 10-15 seconds on an i7/8GB RAM system. With a difference order of 12, it runs for about 20 minutes and consumes about 3-4 gigs of RAM.