Matlab中ode45求解器求解CMUT模型不收敛问题排查

我正在模拟一款可通过带附加非线性力的质量-弹簧-阻尼系统建模的电容微机械超声换能器（CMUT），系统平衡方程如下：
$$\frac{\epsilon_0 a V^2}{2 x^2} - Bx' + Kx = 0$$

我用MATLAB的ode45求解器求解位移x并寻找平衡点，所有系数均按CMUT尺寸定义，代码片段如下：

x0 = 0;
[t, x] = ode45(@(t, x) odefun(t, x, eps0, a, V0, B, d0, K), tspan, x0);

function dxdt = odefun(t, x, eps, a, V, B, d0, K)
    dxdt = ((eps * a * V^2) ./ (2 * B * (d0 - x).^2)) + ((K / B) .* x);
end

方程经导师及多篇文献验证正确，但求解器始终无法收敛，得到的膜位移结果无限增大，这显然不符合CMUT尺寸小于毫米的实际情况。我已尝试MATLAB多种求解器、自行实现Runge-Kutta算法、对照参考文献验证系数，甚至用ode15s刚性求解器并添加额外选项，结果均一致。

请问我使用ode45的方式存在错误吗？

完整代码如下：

% DIMENSIONS DE LA MEMBRANE
e = 500e-9;                 % epaisseur [m]
r = 20e-6;                  % rayon [m]
d0 = 550e-9;                % epaisseur de cavite [m]
a = pi * r^2;               % surface de la membrane [m^2]

% PARAMETRES MECANIQUES
E = 200e9;                  % module d'Young du SiN [Pa]
nu = 0.25;                  % coefficient de poisson du SiN
eta = 18.5e-6;              % viscosité dynamique de l'air [Pa.s]
p0 = 1e5;                   % pression exterieure [Pa]
rhoSiN = 3170;              % densite du SiN [Kg/m^3]
m = rhoSiN * a * e;         % masse membrane [Kg]
K = (16 * E * e^3) / (3 * (1 - nu^2) * r^2);   % raideur [N/m]
B = (eta * pi * r^2) / e;                      % amortissement [N.s/m]

% PARAMETRES ELECTRIQUES
eps0 = 8.85e-12;            % permittivité diélectrique du vide [F/m]
V0 = 10;                    % tension de polarisation [V]

% Conditions initiales et plage de temps
x0 = 0;
ti = 0; tf = 1e-3;
dt = 1e-6;
tspan = linspace(ti, tf, 1/dt);

% Résolution par RK4
[t, x] = ode45(@(t, x) odefun(t, x, eps0, a, V0, B, d0, K), tspan, x0);

plot(t,x,'-')

function dxdt = odefun(t, x, eps, a, V, B, d0, K)
    dxdt = ((eps * a * V^2) ./ (2 * B * (d0 - x).^2)) + ((K / B) .* x);
end

参考文献：

• Y. Wang, L. -M. He, Z. Li, W. Xu and J. Ren, "A computationally efficient nonlinear dynamic model for cMUT based on COMSOL and MATLAB/Simulink"
• T. Merrien, A. Boulmé and D. Certon, "Lumped-Parameter Equivalent Circuit Modeling of CMUT Array Elements"

内容的提问来源于stack exchange，提问作者aloyssss