Having looked around the intergoogles and Mathematica.SE, I thought I'd pose a question with a minimum working example.

Here is the situation I am trying to improve:

1. I am solving a 4th order non linear PDE with NDSolve.
2. It is stiff and I use a stiff solver such as BDF or LSODA.
3. On occassion, I have no choice but to increase the MaxStepFraction to uncomfortable levels.
4. As a result, the code runs longer than usual (made worse by the fact that it is a stiff equation to begin with)

Is there any way I could improve NDSolve performance/speed?

Here is my minimum example:

$HistoryLength = 0;
Needs["VectorAnalysis`"]
Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];
Clear[Eq0, EvapThickFilm, h, Bo, \[Epsilon], K1, \[Delta], Bi, m, r]
Eq0[h_, {Bo_, \[Epsilon]_, K1_, \[Delta]_, Bi_, m_, r_}] := \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]h\) + 
    Div[-h^3 Bo Grad[h] + 
      h^3 Grad[Laplacian[h]] + (\[Delta] h^3)/(Bi h + K1)^3 Grad[h] + 
      m (h/(K1 + Bi h))^2 Grad[h]] + \[Epsilon]/(
    Bi h + K1) + (r) D[D[(h^2/(K1 + Bi h)), x] h^3, x] == 0;
SetCoordinates[Cartesian[x, y, z]];
EvapThickFilm[Bo_, \[Epsilon]_, K1_, \[Delta]_, Bi_, m_, r_] := 
  Eq0[h[x, y, t], {Bo, \[Epsilon], K1, \[Delta], Bi, m, r}];
TraditionalForm[EvapThickFilm[Bo, \[Epsilon], K1, \[Delta], Bi, m, r]];
L = 2*92.389;
 
TMax = 3100*100;
Off[NDSolve::mxsst];
Clear[Kvar];
Kvar[t_] :=  Piecewise[{{1, t <= 1}, {2, t > 1}}]
(*Ktemp = Array[0.001+0.001#^2&,13]*)
hSol = h /. NDSolve[{
      (*Bo,\[Epsilon],K1,\[Delta],Bi,m,r*)
      
      EvapThickFilm[0.003, 0, 1, 0, 1, 0.025, 0],
      h[0, y, t] == h[L, y, t],
      h[x, 0, t] == h[x, L, t],
      (*h[x,y,0] == 1.1+Cos[x] Sin[2y] *)
      
      h[x, y, 0] == 
       1 + (-0.25 Cos[2 \[Pi] x/L] - 0.25 Sin[2 \[Pi] x/L]) Cos[
          2 \[Pi] y/L]
      },
     h,
     {x, 0, L},
     {y, 0, L},
     {t, 0, TMax},
     Method -> {"BDF", "MaxDifferenceOrder" -> 1},
     MaxStepFraction -> 1/50
     ][[1]] // AbsoluteTiming

A BDF limited to Order 1 needs about 41 seconds to solve the equation until failure while the LSODA allowed up to order 12 does a fantastic job of cutting it down to 18 seconds.

Now when I increase the MaxStepFraction, I obviously increase the grid density. I am currently running a case that has several thousand grid points and has taken 24+ HOURS, yes hours and hasn't given me a solution yet. This was expected as I have run cases that took about 3-4 hours before with a coarser grid and do hog the ram (they take up about ~3-4GBs mostly because I am exporting data as .MAT files)

What suggestions could be provided to improve the speed for such a stiff equation?

I have also tried stopping tests and it doesn't quite help all the time as I'd rather mathematica quit my program naturally as a result of overbearing stiffness than artificially through a stopping test. (The former has physical significance)

Yes, this question bears resemblance to this but I don't think its the same.

I have given Parallelize a thought but it doesn't work on NDSolve. Any options that I have either on the Mathematica front, computing front, or saving the interpolation function data?