我正在尝试模拟二元混合物的相分离。如果自由能 F 是浓度的函数，则动力学方程为：$c$

$\frac{\partial c(x,t)}{\partial t}=\frac{d^2}{dx^2} \frac{\delta F[c]}{\delta c}$

对于 Flory-Huggins 自由能，我们有：

$\frac{\delta F[c]}{\delta c}=\log\left(\frac{c}{1-c}\right) + \chi c^2 + \gamma (c')^2$

第一项是熵，第二项是粒子之间的吸引力（），第三项类似于表面张力。$\chi<0$

我使用时间向前和以空间为中心的差异化。即使对于非常小的，我也会得到数值不稳定。我首先认为术语是负责任的，但即使没有，不稳定性仍然存在。$dt$$\gamma(c'^2)$

这是我的无通量边界条件的 Matlab 代码，尽可能清晰和简化。我知道边界条件的实现可能不正确，但我不认为问题来自于此。

我该怎么办？

function phaseSep ()
clear all;
clf;

N = 101;   %number of grid points in x
dx = 1/(N - 1);
x = 0:dx:1;   %vector of x values
T = 1e3;  %number of time steps
dt = 1e-8;

% Second derivative
function y=mylaplace(f,i)
    y = f(i+1)-2*f(i)+f(i-1);
end


function 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 concentration
hold;
plot(x,C);

% iterate
for 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 concentration
plot(x,C);

end