matlab的ode15s解决常微分方程组的问题

问题遇到的现象和发生背景

用matlab的ode15s解一个常微分方程组

问题相关代码，请勿粘贴截图

代码如下（包括一个主函数，两个用于调用的函数Pe、Pi）：
主函数：
function dydt=test(t,y,V,Vp,Vv)
dydt=zeros(5,5);
Pe=pe(y);
Pe=Pe(1:5);
dydt(1)=2/3*Pe./y(4)-y(1)./y(4).dydt(4);%等式右边括号内dydt(4)为离子密度表达式,点运算-->对应元素运算
%电子温度方程
Pi=pi();
%tauE=ta(y);
dydt(2)=2/3Pi./y(4)-y(2)./y(4).*dydt(4);%等式右边括号内dydt(4)为离子密度表达式
%离子温度方程
Rp=rp(y);
L=inductance();
dydt(3)=V-y(3).*Rp/L;
%等离子体电流方程
tauE=ta(y); %约束时间
dydt(4)=y(4)*y(5)*iz-y(4)*y(4)*rec-y(4)tauE;
%氢离子密度
S=speed(y);
dydt(5)=Vp/Vv(y(5)tauE-Sy(5)*y(4));%Lloyd的模型,tauE=tauP
%氢原子密度
end

function Lf=lf(y)%有效连接长度
Bt=2.65;
Bp=2e-3;
a=1.6;
Lf=0.75aBt/Bp*exp(y(3)/100);
end

function tauE=ta(y)%约束时间的倒数
Bt=2.65;
a=1.6;
Cs=cs(y);
Lf=lf(y);
tauE=y(1)/(8a^2Bt)+Cs/Lf;%tauE=1/taue;tanuP=1/taup,约束时间
end

function Cs=cs(y)%离子声速
mi=1.674e-27;
Cs=sqrt((y(1)+y(2))/mi);
end

function Rp=rp(y)%等离子体电阻
R=5.65;
a=1.6;
Rp=(2R/a^2)5.2647e-510y(1)^(-1.5);
end

function L=inductance()
li=0.5;
R=5.65;
a=1.6;
mu=4pi10e-7;
L=muR(log(8*R/a)+li/2-2);
end

function S=speed(y)
%反应速率-Lloyd文章
S=2e-13/(6+y(1)/13.6)*sqrt(y(1)/13.6)*exp(-13.6/y(1));
end

Pe：
function Pe=pe(y)
%电子功率
Poh=poh(y);
Piz=piz(y);
Prad=prad(y);
Pequi=pequi(y);
Pecon=pecon(y);%功率的局部函数，为解决输入参数太多的问题
Pe=Poh-(Piz+Prad)-Pequi-Pecon;
end
function Poh=poh(y)
Vp=200;
Rp=rp(y);
Poh=y(3)^2*Rp/Vp;
end

function Piz=piz(y)
fileID=fopen('E:\test\try\scd96_h.dat');
iz=fread(fileID);%ASAD
iz=iz(1:10097);
energy=13.6;%氢原子的电离能量
Piz=y(4)*y(5)izenergy;
end

function Prad=prad(y)
energy=13.6;
fileID=fopen('E:\test\try\prb96_h.dat') ;
recbr=fread(fileID);
recbr=recbr(1:10097);%令下式矩阵维度一致（1-10097）

fileID=fopen('E:\test\try\plt96_h.dat');
line=fread(fileID);
line=line(1:10097);%矩阵维度

fileID=fopen('E:\test\try\acd96_h.dat');
rec=fread(fileID);%ADSA
rec=rec(1:10097);%令维度一致

Prad=y(5)*y(4)*line+y(4)*y(4)*recbr-y(4)*y(4)recenergy;
end

function Pequi=pequi(y)
Pequi=7.75e-34*y(4)*y(4)10(y(1)-y(2))/y(1)^1.5;
end

function Pecon=pecon(y)
tauE=ta(y);
Pecon=2/3*y(4)*y(1)*tauE;
end

function Rp=rp(y)%等离子体电阻
R=5.65;
a=1.6;
Rp=(2R/a^2)5.2647e-510y(1)^(-1.5);
end

function tauE=ta(y)%约束时间的倒数
Bt=2.65;
a=1.6;
Cs=cs(y);
Lf=lf(y);
tauE=y(1)/(8a^2Bt)+Cs/Lf;%tauE=1/taue;tanuP=1/taup,约束时间
end

function Cs=cs(y)%离子声速
mi=1.674e-27;
Cs=sqrt((y(1)+y(2))/mi);
end

function Lf=lf(y)%有效连接长度
Bt=2.65;
Bp=2e-3;
a=1.6;
Lf=0.75aBt/Bp*exp(y(3)/100);
end

Pi：

function Pi=pi(y)
%离子功率
Pequi=pequi(y);
Pcx=pcx(y);
Pconv=pconv(y);
Pi=Pequi-Pcx-Pconv;
end

function Pequi=pequi(y)
Pequi=7.75e-34*y(4)*y(4)10(y(1)-y(2))/y(1)^1.5;%10为库伦对数
end

function Pcx=pcx(y)
fileID=fopen('E:\test\try\ccd96_d.dat');
cx=fread(fileID);
T0=0.026;
n0=4.8e200.8e-3;%p0=0.8mPa
Pcx=1.5(y(2)-T0)n0y(4)*cx;
end

function Pconv=pconv(y)
tauE=ta(y);
Pconv=3/2*y(4)*y(2)*tauE;
end

function tauE=ta(y)%约束时间的倒数
Bt=2.65;
a=1.6;
Cs=cs(y);
Lf=lf(y);
tauE=y(1)/(8a^2Bt)+Cs/Lf;%tauE=1/taue;tanuP=1/taup,约束时间
end

function Cs=cs(y)%离子声速
mi=1.674e-27;
Cs=sqrt((y(1)+y(2))/mi);
end

运行函数的文件：
V=12;
Vp=200;
Vv=1000;
tspan=[0 0.3];
Te0=1;%y(1)
Ti0=0.026;%y(2)
I0=2.4e3;%y(3)
nH1=0;%y(4)
p0=0.8e-3;%气压
nH0=4.8e20*p0;%y(5)
y0=[Te0,Ti0,I0,nH1,nH0];
[t,y]=ode15s(@test,tspan,y0);
plot(t,y);

运行结果及报错内容

我的解答思路和尝试过的方法

最初，等式两边的类型、维度都查看了，发现都是double，只是维度不同，我将乘除用点乘和点除替换。可还是出现了上述问题。现在，我认为是dydt=zeros(5,1);有问题。

我想要达到的结果

使ode15s能够运行，函数不出错。