# Numerical instability when using NDSolve to simulate phase separation

I'm using NDSolve to simulate the following model of phase separation:

where rho_T=rho_1+rho_2. The term with chi causes rho_1 and rho_2 to repel each other. However, there is a numerical instability at the interface which causes the solution to blow up:

Refining the mesh narrows the cusps, but they still diverge and cause the integration to halt. How can I use NDSolve to solve this system?

Here is the code in Mathematica 13.1. Note that decreasing chi allows NDSolve to converge, but you end up with a uniform solution instead of two domains.

rhoT = rho1[t, x] + rho2[t, x];

chi = 10.5;

vars = {{rho1[t, x], rho2[t, x]}, t, {x}};
op = D[vars[[1]], t] +
DiffusionPDETerm[
vars, {{alpha1/rhoT,
alpha1*(rho1[t, x]/rhoT)*chi}, {alpha2*(rho2[t, x]/rhoT)*chi,
alpha2/rhoT}}];

(*no-flux boundary conditions*)

bcConserved = NeumannValue[0, x == 0] + NeumannValue[0, x == l];

alpha1 = 1;
alpha2 = 1;

tF = 1;
e0 = 0.1;
l = 1;

(*initial conditions*)

f1[x_] := (e0/
l)*(20/(Log[20*(1 + E^(10*l))] - Log[20*(1 + E^(-10*l))]))/(1 +
E^(20*(x - l/2)));

f2[x_] := (e0/
l)*(-(20/(Log[20*E^(10*l) (1 + E^(10*l))] -
Log[20*(1 + E^(10*l))]))/(1 + E^(20*(x - l/2))) + 2);

ic1 = rho1[0, x] == f1[x];

ic2 = rho2[0, x] == f2[x];

sol = NDSolveValue[{op == {bcConserved, bcConserved}, ic1,
ic2}, {rho1, rho2}, {t, 0, tF}, {x, 0, l}];

{r1, r2} = sol;

Plot[{r1[tF, x], r2[tF, x]}, {x, 0, l}, PlotLegends -> {"1", "2"},
PlotRange -> {Automatic, Full}, LabelStyle -> Directive[Black, 10]]


• This is a typical problem with pure Neumann boundary conditions discussed on mathematica.stackexchange.com/questions/191476/… Jan 13 at 4:33
• @AlexTrounev, how do get to this conclusion? This is a time dependent problem, while the one you link to is a stationary problem. Jan 13 at 5:52
• The second equation is independent from the first, is that right? Where does this equation come from? Jan 13 at 5:57
• @user21 There is a typo in the Latex form equation for $\rho_2$. Equations should be symmetrical. But in the code there is the symmetry. Jan 13 at 6:48
• @user21 It's a mass-transfer equation where rho is the mass density and the flux J=-alpha * concentration * gradient of chemical potential. The concentration is rho/rhoT. The chemical potential is mu_i =ln(rho_i)+chi * rho_j. (I'm playing fast and loose with units here, but that's not the source of the problem.) It's very similar to Eq. 10 in this paper, although I'm working on a slightly different problem. pubs.rsc.org/en/content/articlelanding/2019/sm/c8sm02045k Jan 13 at 16:37