背景
我有一个 6 个 ODE 的刚性系统，在 MATLAB 中表示如下：
system = @(t,x)[x(1,:).*x(2,:).*(-5.726882618492327e8)-x(1,:).*x(3,:).*1.467710449114545e10-x(1,:).*x(4,:).*3.012162507288153e12+x(3,:).*x(6,:).*1.137674873581712e12+x(4,:).*x(6,:).*5.703626484852603e12+x(5,:).*x(6,:).*1.656793950858653e9-x(1,:).^2.*2.980001443206009e8+x(6,:).^2.*8.25e9;x(1,:).*x(2,:).*(-5.726882618492327e8)+x(3,:).*x(6,:).*1.137674873581712e12-x(2,:).^2.*4.476510081133096e2+x(6,:).^2.*3.403491297089207e3;x(1,:).*x(2,:).*5.726882618492327e8-x(1,:).*x(3,:).*1.467710449114545e10-x(3,:).*x(6,:).*1.137674873581712e12+x(4,:).*x(6,:).*5.703626484852603e12;x(1,:).*x(3,:).*1.467710449114545e10-x(1,:).*x(4,:).*3.012162507288153e12-x(4,:).*x(6,:).*5.703626484852603e12+x(5,:).*x(6,:).*1.656793950858653e9;x(5,:).*(-8.186979330234089e8)+x(6,:).*2.1e9+x(1,:).*x(4,:).*3.012162507288153e12-x(5,:).*x(6,:).*1.656793950858653e9;x(5,:).*8.186979330234089e8-x(6,:).*2.1e9+x(1,:).*x(2,:).*5.726882618492327e8+x(1,:).*x(3,:).*1.467710449114545e10+x(1,:).*x(4,:).*3.012162507288153e12-x(3,:).*x(6,:).*1.137674873581712e12-x(4,:).*x(6,:).*5.703626484852603e12-x(5,:).*x(6,:).*1.656793950858653e9+x(1,:).^2.*2.980001443206009e8+x(2,:).^2.*4.476510081133096e2-x(6,:).^2.*8.250003403491297e9];
如果我使用 ODE 求解器，我可以求解 ODE（此处显示）并使用它来近似稳态解。举个例子，
[t,x] = ode15s(system,[0 1e-6],[0 0 0 0 0 1]);
plot(t,x)
ss_sol = x(end,:);
作为参考，这样得到的稳态解近似为
ss_sol = [0.322330352943109, 0.458043435766086, 0.001213186698607, 0.000016426443105, 0.157136142002196, 0.061260531336451];
即使我稍微改变初始条件，我也会回到这个解决方案。请注意，sum(ss_sol)大约为 1（ODE 是在此假设下得出的）。另请注意，由于物理限制，每个变量必须介于 0 和 1 之间，解决方案满足此条件。
问题
为了获得稳态解，我想将导数设置为零，并为其根求解类似的非线性代数系统。该系统在 MATLAB 中稍作改写为
system = @(x)[x(1).*x(2).*(-5.726882618492327e8)-x(1).*x(3).*1.467710449114545e10-x(1).*x(4).*3.012162507288153e12+x(3).*x(6).*1.137674873581712e12+x(4).*x(6).*5.703626484852603e12+x(5).*x(6).*1.656793950858653e9-x(1).^2.*2.980001443206009e8+x(6).^2.*8.25e9,x(1).*x(2).*(-5.726882618492327e8)+x(3).*x(6).*1.137674873581712e12-x(2).^2.*4.476510081133096e2+x(6).^2.*3.403491297089207e3,x(1).*x(2).*5.726882618492327e8-x(1).*x(3).*1.467710449114545e10-x(3).*x(6).*1.137674873581712e12+x(4).*x(6).*5.703626484852603e12,x(1).*x(3).*1.467710449114545e10-x(1).*x(4).*3.012162507288153e12-x(4).*x(6).*5.703626484852603e12+x(5).*x(6).*1.656793950858653e9,x(5).*(-8.186979330234089e8)+x(6).*2.1e9+x(1).*x(4).*3.012162507288153e12-x(5).*x(6).*1.656793950858653e9,x(5).*8.186979330234089e8-x(6).*2.1e9+x(1).*x(2).*5.726882618492327e8+x(1).*x(3).*1.467710449114545e10+x(1).*x(4).*3.012162507288153e12-x(3).*x(6).*1.137674873581712e12-x(4).*x(6).*5.703626484852603e12-x(5).*x(6).*1.656793950858653e9+x(1).^2.*2.980001443206009e8+x(2).^2.*4.476510081133096e2-x(6).^2.*8.250003403491297e9];
但是，我很难获得ss_sol前面提到的相同值。
使用vpasolve，即使我将精度提高到 100+ 位，也不会导致任何解决方案。使用fsolveor lsqnonlin（后者让我分别设置 0 和 1 的下限和上限）收敛到不是稳态解的值，除非我提供ss_sol作为初始猜测。我尝试降低公差并增加迭代次数，但仍然没有运气。
有什么建议？将 ODE 集成到稳态的最佳方法真的是？我一直在阅读，将 ODE 的刚性系统简化为代数方程应该更容易求解，但我无法以这种方式找到稳态解。我需要一种通用方法来解决许多不同系统的稳态值，因此问题超出了上述具体示例。