function dydz=TM1CR(z,y)
global N1 N0
N=5.51E24;%掺铥光纤Tm浓度
t1=334.7E-6;
t3=14.2E-6;
c=3E8;
beta30=0.14;
beta31=0.72;
A10=1/t1;
A30=beta30/t3;
A31=beta31/t3;
T1=0.8964;
L=1;
h=6.63E-34;
A=1.39E-11;
R1=0.05;
R2=0.04;
R3=1;
R4=0.04;
darp=5E-25;
dars=10E-26;
ders=2.5E-25;
k1=1.8E-22;
k2=0.084*k1;
dp=1.2E-2;
ds=2.3E-3;

rse=c*(y(3)+y(4))*ders;
rsa=c*(y(3)+y(4))*dars;
rp=darp*(y(1)+y(2));
A1=A10+T1+rse-2*A30-A31;
A2=2*rp+rsa+2*A30+A31;
N1=-1/2*(-k1*N*A2*A1+2*k1*A2*N*A30+k1*A2*N*A31+2*k1*N*A31*A1+4*k1*N*A30*A1-(k1^2*N^2*A2^2*A1^2+4*k1^2*A2^2*N^2*A30^2+k1^2*A2^2*N^2*A31^2+4*k1^2*N^2*A2^2*A1*A30+2*k1^2*N^2*A2^2*A1*A31+4*k1^2*A2^2*N^2*A30*A31-4*k1*A2^3*A1*k2*N1^2-4*k1*A1^2*k2*N1^2*A2^2)^(1/2))/k1/A1/(A2+A1);
N0=((2*A30+A31)*N-1/2*(-k1*N*A2*A1+2*k1*A2*N*A30+k1*A2*N*A31+2*k1*N*A31*A1+4*k1*N*A30*A1-(k1^2*N^2*A2^2*A1^2+4*k1^2*A2^2*N^2*A30^2+k1^2*A2^2*N^2*A31^2+4*k1^2*N^2*A2^2*A1*A30+2*k1^2*N^2*A2^2*A1*A31+4*k1^2*A2^2*N^2*A30*A31-4*k1*A2^3*A1*k2*N1^2-4*k1*A1^2*k2*N1^2*A2^2)^(1/2))/k1/(A2+A1))/A2;

dydz=[-y(1)*(darp*N0+dp)
    y(2)*(darp*N0+dp)
    y(3)*(ders*N1-dars*N0-ds)
    -y(4)*(ders*N1-dars*N0-ds)];
以上存为TM1CR.m文件
plaunched=5;
R1=0.05;R2=0.04;R3=1;R4=0.04;
pr0=0.01;sr0=1.5;
pf0=plaunched+R1*pr0;
% prL=R2*pfL;
sf0=R3*sr0;
% srL=R4*sfL
x0=[pf0;pr0;sf0;sr0];
[z,y]=ode45(@TM1CR,[0,1],x0);
subplot(2,2,1);
plot(z,y(:,1),'-.',z,y(:,2),'-',z,y(:,1)+y(:,2),'*');

此部分调用上面的.m文件，但是一直报错，如下：
??? Error using ==> odearguments at 117
Solving TM1CR requires an initial condition vector of length 4.

Error in ==> ode45 at 173
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Error in ==> pfprcom at 10
[z,y]=ode45(@TM1CR,[0,1],x0);

针对第一个错误提示，我觉得我的初始条件是长度为4的向量，可是却一直报错，一直调不好，求达人帮助，感激不尽！多谢！