MATLAB实现文章化学方程组求解结果不符，求问题排查

我最近在复现论文《The chemistry of co-injected BOE》里的化学物种占比计算，把论文里的化学方程组写成了MATLAB代码，但求解结果和论文图2的参考值一直对不上，反复核对数据也没找到问题，真心求大佬帮忙排查一下！

论文里的目标参考图是图2：

我在10cc、40cc、90cc三种条件下，求解结果和论文参考值的对比情况如下：

我用到的MATLAB主脚本代码如下：

clc; clear all; 
%Units to be used 
%Volume is in CC also cm^3, 1 litre is 1000 CC, 1 cc = 1 ml 
%density is in g/cm^3 
%weigth percentages are in fractions of 0 to 1 
%Molecular weight is in g/mol 

% pts=10; %number of points for linear spacing 
%weight percentages of NH4OH and HF 
xhf=0.49; 
xnh3=0.28; 

%H2O 
Vh2o=1800; 
dh2o=1.00; %0.997 at 25C when rounded 1 
mh2o=18.02; 

%HF values 
Vhf=100; 
dhf49=1.15; 
dhf=dh2o+(dhf49-dh2o)*xhf/0.49; %@ 25C 
Mhf=20.01; 
nhf=mols(Vhf,dhf,xhf,Mhf); 

%NH4OH (NH3) values 
% Vnh3=linspace(0.1*Vhf,1.9*Vhf,pts); 
Vnh3=10; 
dnh3=0.9; %for ~20-31% @~20-25C 
Mnh3=17.03; %The wt% of NH4OH actually refers to the wt% of NH3 dissolved in H2O 
nnh3=mols(Vnh3,dnh3,xnh3,Mnh3); 

if max(nnh3)>=nhf 
    error(['There are more mols NH4OH,',num2str(max(nnh3)),', than mols HF,',num2str(nhf),'.']) 
end 

%% Calculations for species 
Vt=(Vhf+Vh2o+Vnh3)/1000; %litre 
A=nhf/Vt; %mol/l 
B=nnh3/Vt; %mol/l 

syms HF F H2F2 HF2 NH3 NH4 H OH 

eq2= H*F/HF==6.85*10^(-4); 
eq3= NH3*H/NH4==6.31*10^(-10); 
eq4= H*OH==10^(-14); 
eq5= HF2/(HF*F)==3.963; 
eq6= H2F2/(HF^2)==2.7; 
eq7= H+NH4==OH+F+HF2; 
eq8= HF+F+2*H2F2+2*HF2==A; 
eq9= NH3+NH4==B; 

eqns=[eq2,eq5,eq6,eq8,eq4,eq3,eq9,eq7]; 
varias=[HF, F, H2F2, HF2, NH3, NH4, H, OH]; 

assume(HF> 0 & F>= 0 & H2F2>= 0 & HF2>= 0& NH3>= 0 & NH4>= 0 & H>= 0 & OH>= 0) 

[HF, F, H2F2, HF2, NH3, NH4, H, OH]=vpasolve(eqns,varias);% [0 max([A,B])]) 

totalHF=double(HF)+double(F)+double(H2F2)+double(HF2); 
HFf=double(HF)/totalHF %fraction of species for HF 
H2F2f=double(H2F2)/totalHF %fraction of species for H2F2 
Ff=double(F)/totalHF %fraction of species for F- 
HF2f=double(HF2)/totalHF %fraction of species for HF2- 

还需要一个额外的mols.m辅助函数：

%%%% amount of mol, Vol=volume, d=density, pwt=%weight, M=molecularweight 
function mol=mols(Vol, d, pwt, M) 
    mol=(Vol*d*pwt)/M; 
end 

论文里用到的化学平衡方程式如下：

注：我的脚本里的H2F2对应论文中的(HF)₂。

内容的提问来源于stack exchange，提问作者Bob van de Voort