matlab中用ode23画图，为什么函数一直报错

matlab中用ode23画图，为什么函数一直报错

调用的函数

function dy=odefun(t,y)
global s d beta p  k1 c elta k2 m1 m alpha delta b
dy=zeros(4,1);    
dy(1)=s-d.*y(1)-beta.*y(1).*y(2);
dy(2)=k1.*y(1).*y(2)+alpha.*y(3)-delta.*y(2)-p.*y(4).*y(2);
 dy(3)=k2.*y(1).*y(2)-m1.*y(3);
 dy(4)=((c.*y(2).*y(4))./(1+elta.*y(2)))-m.*y(2).*y(4)-b.*y(4);

程序


clear %%Et后向分支
clc
%global s d beta p c elta k1 k2 m1 m alpha delta 

 s=10000;   d=0.01;  p=1.8*10^(-8); beta=2.2*10^(-4); elta=0.01;  b=60;  c=20;  m=8; k1=2.4*10^(-8); delta=0.01; k2=2.4*10^(-8);alpha=1.7946;m1=1.2946;

% alpha=(-4)./5.*(exp(-0.05.*100)-1)+0.01.*100;
%m1=(-4)./5.*(exp(-0.05.*100)-1)+0.5.*100*0.01;
 v1=m+b.*elta-c;
 u1=4.*b.*m.*elta;
 Delta1=v1^2-u1;
c1=roundn(b.*elta+m+2.*sqrt(m.*elta.*b),-10); 
 %alpha=8.8*((-exp(-10000)+0.5.*exp(0.25-20000))-(-exp(-0.25)+0.5.*exp(0.25-0.5)));
 Q=m1./alpha;
 R0=(s.*(k1+k2.*Q))./(d.*delta);
 T1=s./(d.*R0);
 I1=(d.*(R0-1))./beta;
 E1=0;
I2=sqrt(b./(m.*elta));
T2=s./(d+beta.*I2);
E2=((k1+k2.*Q).*T2-delta)./p;
 L1=k2.*T1.*I1;

c2=roundn(m+b.*elta+m.*elta.*d.*(R0-1)./beta+b.*beta./(d.*(R0-1)),-10);
 I_=(c-m-(b.*elta)-sqrt((b.*elta+m-c)^2-(4.*m.*elta.*b)))./(2.*m.*elta);%I星负
 I__=(c-m-b.*elta+sqrt((b.*elta+m-c)^2-4.*m.*elta.*b))./(2.*m.*elta);%I星正
 T_=s./(d+beta.*I_);
 E_=(T_.*(k1+k2.*Q)-delta)./p;%E星负

 T__=s./(d+beta.*I__);
 E__=(T__.*(k1+k2.*Q)-delta)./p;%E星正
% v=T__.*(k1+k2.*Q)-delta;
 %v3=T__.*(k1+k2.*Q);
 L2=k2.*T_.*I_;
 Rc=1+(beta.*sqrt(b./(m.*elta)))./d;
 R21=s.*(k1+k2.*Q)./(delta.*(d+beta.*I_));%R星负
 R11=s.*(k1+k2.*Q)./(delta.*(d+beta.*I__));%R星正
% A=s.*(k1+k2.*Q);
% B=delta.*(d+beta.*I__);

options=odeset('RelTol',1e-10,'AbsTol',[1e-10 1e-10 1e-10 1e-10]);
%[t,y]=ode23('odefun', [0 600],[10^6 0 0 0.001]);
[t,y]=ode23('odefun',[0 500],[10^6 0 0 0.001],options);
% E1=（3，0） 
% E+=(2,1)
% E-=(0.5，2)

figure(1);
plot(t,y(:,1),'r','LineWidth',2);xlabel('Days');ylabel('T/T cells')
set(gca,'FontSize',20)
hold on

figure(2) ;
plot(t,y(:,2),'r','LineWidth',2);xlabel('Days');ylabel('I/I virus')
set(gca,'FontSize',20)
hold on
% 
 figure(3) ;
 plot(t,y(:,3),'r','LineWidth',2);xlabel('Days');ylabel('L/L(t,a)')
 set(gca,'FontSize',20)
 hold on
 
 figure(4) ;
  plot(t,y(:,4),'r','LineWidth',2);xlabel('Days');ylabel('E/Immune cells')
 set(gca,'FontSize',20)
 hold on

结果一直报错，如下图：

img