I am trying to simulate the phase separation of a binary mixture. If the free energy F is known as a function of the concentration $c$, the dynamical equation is:$\frac{\partial c(x,t)}{\partial t}=\frac{d^2}{dx^2} \frac{\delta F[c]}{\delta c}$For the Flory-Huggins free energy we have:$\frac{\delta F[c]}{\delta c}=\log\left(\frac{c}{1-c}\right) + \chi c^2 + \gamma (c')^2$First term is entropy, second term is attraction between particles ($\chi<0$) and third term is similar as a surface tension.I use time forward and space centered differentiation. Even for very small $dt$ I get numerical instabilities. I first thought the term $\gamma(c'^2)$ was responsible, but instabilities remain even without.Here is my Matlab code for no-flux boundary conditions, clear and simplified as much as I can. I am aware the boundary condition implementation may not be correct but I don't think the problem comes from this.What should I do?function phaseSep ()clear all;clf;N = 101;   %number of grid points in xdx = 1/(N - 1);x = 0:dx:1;   %vector of x valuesT = 1e3;  %number of time stepsdt = 1e-8;% Second derivativefunction y=mylaplace(f,i)    y = f(i+1)-2*f(i)+f(i-1);endfunction y=dfdc(C)    p=-0.01;    g=0.01;    for i=2:N-1    y(i) = log(C(i)/(1-C(i)))+p*C(i)+g*mylaplace(C,i)/dx^2;    end    % No-flux Boundary conditions    y(1)=y(2);    y(N)=y(N-1);end%Initial concentration  for i=1:N    C(i)=0.2+0.1*tanh(10*(x(i)-0.5));end% Plot initial concentrationhold;plot(x,C);% iteratefor t=1:T    Z=dfdc(C);   for i=2:N-1      C(i) = C(i) + (mylaplace(Z,i)/dx^2)*dt;    end   %No flux boundary conditions   C(1) = C(2);   C(N) = C(N-1);end% Plot final concentrationplot(x,C);end
Modified diffusion equation and unstabilities
finite difference;diffusion
At first glance, it looks like you are using the method of lines with forward Euler time steps. What WolfgangBangerth is getting at with his example is that even for a simple heat equation, the stability limit of forward Euler (namely, that $|\lambda\Delta{t}| < 1$) combined with the eigenvalues induced by a finite difference approximation of the diffusion operator result in a very stringent limit on the time step relative to the spatial discretization. For a second-order centered finite difference approximation to the second-order derivative term, this limit is of the form $\Delta{t} < h^{2}/(2a)$, where $a$ is the diffusivity and $h$ is the grid spacing of the spatial discretization (assumed uniform). You should check out a text such as LeVeque's Finite Difference Methods for Ordinary and Partial Differential Equations for a basic overview of the theory. In terms of practical advice, you probably want to replace the time discretization you are currently using with an L-stable implicit method. To start, you might try using ode23tb in MATLAB.