Solving simultaneous ODEs with NDSolve causes unexpected singularity error (ndsz)

I’m trying to model the evaporation of a binary droplet in a square well by solving a couple of PDEs for its height profile and composition. I can manage it for a single liquid pretty well, but when I introduce the PDE for composition NDSolve quickly finds a singularity in the composition variable (importantly the singularity moves slightly if I change the number of points). I can’t seem to find any way to stabilise the system, can anyone help me? Is there something fundamental about NDSolve that I’m missing?

(I asked a similar question a while ago but I’ve restructured my code to make it easier to share here)

The equations are:

$ $ \frac{\partial h}{\partial t} = – C \frac{\partial}{\partial x}\Big[\frac{1}{3} \frac{\partial^3h}{\partial x^3}h^3 + \frac{B1}{2} \frac{\partial X}{\partial x}\Big] – (1-AX) $ $ $ $ \frac{\partial X}{\partial t} = -\frac{C}{3} h^2 \frac{\partial^3h}{\partial x^3}\frac{\partial X}{\partial x} – \frac{hB1}{2} \Big(\frac{\partial X}{\partial x}\Big)^2 + \frac{A}{h}X(1-X) $ $

with height $ h$ and composition $ X$ .

And the boundary conditions are

$ h(t)=1$ and $ X(t)=X(t=0)$ at $ x=1$ , $ \frac{\partial h}{\partial x}=\frac{\partial X}{\partial x}=0$ at $ x=0$ , and $ \frac{\partial^3h}{\partial x^3}=0$ at $ x=0$ and $ x=1$ .

I discretise all this using finite differences.

This is what I get for a single liquid droplet (when I suppress equation 2) and I expect something qualitatively similar for a binary droplet. Droplet curves when the composition equation is suppressed

Here is a minimum replicable example:

    N1 = 100; dx = 1/N1; (*Discretise*)     C1  = 1; (*parameter C*)     a = 1.5; b = 1 - a; hInitial = Table[a + b i^2 dx^2, {i, 0, N1}]; (*Initial height profile*)     E1 = 1;(*evaporation on/off switch*)     Cc = C1/(24 * dx^3);     A1 = 0.01; B1 = 1;  (*Control constants in the equations above*)     CcX = (C1 B1)/(2 dx^2); CcX2 = C1/(6 dx^2);     X0 = 0.4; XInitial = Table[X0, {i, 0, N1}]; (*Initial composition ratio*)     Sc = 1; (*Composition on/off switch*)     TGap = 0.085;(*for spacing curves in the plot*)         dv[v_List] :=       With[{h = Take[v, Length[v]/2 - 1], hN = v[[(Length[v]/2)]],          X1 = v[[Length[v]/2 + 1 ;; Length[v] - 1]], XN = v[[-1]]},        With[{          dh1 =            ListCorrelate[{0, 0, 1 , 1, 0}, #] &, (*derivatives are discretised as finite differences.            ListCorrelate achieves this.*)           dh2 = ListCorrelate[{0, 1, 1, 0, 0}, #] &,          dh3 = ListCorrelate[{0, -1, 3, -3, 1}, #] &,          dh4 = ListCorrelate[{-1, 3, -3, 1, 0}, #] &,          hi = ListCorrelate[{0, 0, 1, 0, 0}, #] &,          dh5 = ListCorrelate[{-2, 1, 0, -1, 2}, #] &,            dx1 = ListCorrelate[{0, -1, 1}, #] &,          dx2 = ListCorrelate[{-1, 1, 0}, #] &,          Xi = ListCorrelate[{0, 1, 0}, #] &          },         Flatten@{             -Cc*(dh1[#1]^3*dh3[#1] - dh2[#1]^3*dh4[#1]) -               CcX*(dh1[#1]^2*dx1[#2] - dh2[#1]^2*dx2[#2]) -              E1*(1 - A1 Xi[#2]), (*ODE for heights*)              0,(*height is pinned at edge, (1, 1)*)              Sc*(-CcX2* hi[#1]^2*dh5[#1]*dx2[#2] -              CcX/C1* hi[#1]*dx2[#2]^2 +              A1 *(hi[#1])^-1*Xi[#2]*(1 - Xi[#2])), (*ODE for composition*)                      0 (*composition ratio is fixed at edge, XN = 0.4*)             } &[          Flatten@Join[{h[[3]]}, {h[[2]]}, {h}, {1}, {3 - 3 h[[-1]] + h[[-2]]}],                       Flatten@Join[{X1[[2]]}, {X1}, {XN}]]]                (*this contains the boundary conditions*)         ];      v0 = Join[hInitial, XInitial]; (*Initial conditions*)     system2d =       NDSolve[{v'[t] == dv[v[t]], v[0] ==  v0,          WhenEvent[          Min@Table[             Take[values[[tt]], Length[values[[tt]]]/2],               {tt, 1,  Length[values]}] < 0,           "StopIntegration"](*WhenEvent to stop integration at touchdown*)},          v, {t, 0, 2}][[1, 1, 2]];       Needs@"DifferentialEquations`InterpolatingFunctionAnatomy`";     values = InterpolatingFunctionValuesOnGrid@system2d;      times = Flatten@InterpolatingFunctionGrid@system2d;      T1 = 0; timesteps = {}; For[i = 1, i < Length[times], i++,       If[times[[i]] > T1, {T1 = T1 + TGap,         timesteps = Append[timesteps, {i, times[[i]]}]}]];     finaltime = {{-1}}; firstandfinaltime = {{1}, {-1}};      heights =        Table[Take[values[[tt]], Length[values[[tt]]]/2],          {tt, 1, Length[values]}];      comps = Table[        Take[values[[tt]], -(Length[values[[tt]]]/2)],           {tt, 1, Length[values]}];      Show[Table[       ListPlot[Table[{i*dx, heights[[ttt[[1]], i + 1]]}, {i, 0, N1}],         PlotRange -> {{0, 1}, {0, 1.05*a}}, Joined -> True,         AxesLabel -> {x, h}], {ttt, finaltime}]]      Show[Table[       ListPlot[Table[{i*dx, comps[[ttt[[1]], i + 1]]}, {i, 0, N1}],         PlotRange -> {{0, 1}, {0, 1}}, Joined -> True,        AxesLabel ->{x, X}], {ttt, finaltime}]] 

Thank you for any help you can offer me!