function y=ZYKZ(s1,p,a,c1,d1,a1,vc,r,b,c3,a2,d2,x0,y0,z0) %建立关于模型中参数的函数
%ZYKZ(1.3*10^4,0.1254,2.02*10^7,3.42*10^(-10),0.0412,0.034,0.063,0.18,2*10^(-9),0.006,0.9,0.3466,1,5*10^6,1)
test=-1;%收敛测试变量，17行运行while,一旦发生收敛则变为非负
delta=0.001; N=10;
t=linspace(0,100,N+1);%时间长度(0，1)，等距节点有1001个，包括0，1
h=1/N;%间距
h2=h/2;%用于龙格库塔子程序
u=zeros(1,N+1);
v=zeros(1,N+1);
x=zeros(1,N+1);
x(1)=x0;
y=zeros(1,N+1);
y(1)=y0;
z=zeros(1,N+1);
z(1)=z0;
lambda1=zeros(1,N+1);%建立三个变量向量，每个时间点对应值，均为1行1001列
lambda2=zeros(1,N+1);
lambda3=zeros(1,N+1);
while(test<0)
oldu=u;
oldv=v;
oldx=x;
oldy=y;
oldz=z;
oldlambda1=lambda1;%收敛测试需要当前和以前迭代的值，将当前的值诸存为以前的值
oldlambda2=lambda2;
oldlambda3=lambda3;
for i=1:N
k1=s1+p*x(i)*y(i)/(a+y(i))-c1*x(i)*y(i)-d1*x(i)-a1*(1-exp(-vc*z(i)))*x(i)+u(i); %x状态方程的右端
k2=s1+p*(x(i)+h2*k1)*(y(i)+h2*k1)/(a+(y(i)+h2*k1))-c1*(x(i)+h2*k1)*(y(i)+h2*k1)-d1*(x(i)+h2*k1)-a1*(1-exp(-vc*(z(i)+h2*k1)))*(x(i)+h2*k1)+0.5*(u(i)+u(i+1));
k3=s1+p*(x(i)+h2*k2)*(y(i)+h2*k2)/(a+(y(i)+h2*k2))-c1*(x(i)+h2*k2)*(y(i)+h2*k2)-d1*(x(i)+h2*k2)-a1*(1-exp(-vc*(z(i)+h2*k2)))*(x(i)+h2*k2)+0.5*(u(i)+u(i+1));
k4=s1+p*(x(i)+h*k3)*(y(i)+h*k3)/(a+(y(i)+h*k3))-c1*(x(i)+h*k3)*(y(i)+h*k3)-d1*(x(i)+h*k3)-a1*(1-exp(-vc*(z(i)+h*k3)))*(x(i)+h*k3)+u(i+1);
x(i+1)=x(i)+(h/6)*(k1+2*k2+2*k3+k4);
end
for i=1:N
k1=r*y(i)*(1-b*y(i))-c3*x(i)*y(i)-a2*(1-exp(-vc*z(i)))*y(i);%y方程状态方程的右端
k2=r*(y(i)+h2*k1)*(1-b*(y(i)+h2*k1))-c3*(x(i)+h2*k1)*(y(i)+h2*k1)-a2*(1-exp(-vc*(z(i)+h2*k1)))*(y(i)+h2*k1);
k3=r*(y(i)+h2*k2)*(1-b*(y(i)+h2*k2))-c3*(x(i)+h2*k2)*(y(i)+h2*k2)-a2*(1-exp(-vc*(z(i)+h2*k2)))*(y(i)+h2*k2);
k4=r*(y(i)+h*k3)*(1-b*(y(i)+h*k3))-c3*(x(i)+h*k3)*(y(i)+h*k3)-a2*(1-exp(-vc*(z(i)+h*k3)))*(y(i)+h*k3);
y(i+1)=y(i)+(h/6)*(k1+2*k2+2*k3+k4);
end
for i=1:N
k1=-d2*z(i)+v(i);%z方程状态方程的右端
k2=-d2*(z(i)+h2*k1)+0.5*(v(i)+v(i+1));
k3=-d2*(z(i)+h2*k2)+0.5*(v(i)+v(i+1));
k4=-d2*(z(i)+h*k3)+v(i+1);
z(i+1)=z(i)+(h/6)*(k1+2*k2+2*k3+k4);
end
for i=1:N
j=N+2-i;
k1=lambda1(j)*(c1*y(j)+d1+a1*(1-exp(-vc*z(j)))-p*y(j)/(a*y(j)))+lambda2(j)*c3*y(j);%伴随变量的微分方程，伴随公式
k2=(lambda1(j)-h2*k1)*(c1*0.5*(y(j)+y(j-1))+d1+a1*(1-exp(-vc*0.5*(z(j)+z(j-1))))-p*0.5*(y(j)+y(j-1))/(a*0.5*(y(j)+y(j-1))))+(lambda2(j)-h2*k1)*c3*0.5*(y(j)+y(j-1));
k3=(lambda1(j)-h2*k2)*(c1*0.5*(y(j)+y(j-1))+d1+a1*(1-exp(-vc*0.5*(z(j)+z(j-1))))-p*0.5*(y(j)+y(j-1))/(a*0.5*(y(j)+y(j-1))))+(lambda2(j)-h2*k2)*c3*0.5*(y(j)+y(j-1));
k4=(lambda1(j)-h*k3)*(c1*0.5*y(j-1)+d1+a1*(1-exp(-vc*0.5*z(j-1)))-p*0.5*y(j-1)/(a*0.5*y(j-1)))+(lambda2(j)-h*k3)*c3*0.5*y(j-1);
lambda1(j-1)=lambda1(j)-(h/6)*(k1+2*k2+2*k3+k4);
end
for i=1:N
j=N+2-i;
k1=lambda1(j)*(c1*x(j)-p*a*x(j)/(a+y(j))^2)+lambda2(j)*(2*r*b*y(j)+c3*x(j)+a2*(1-exp(-vc*z(j)))-r);%伴随变量的微分方程，伴随公式
k2=(lambda1(j)-h2*k1)*(c1*0.5*(x(j)+x(j-1))-p*a*0.5*(x(j)+x(j-1))/(a+0.5*(y(j)+y(j-1)))^2)+(lambda2(j)-h2*k1)*(2*r*b*0.5*(y(j)+y(j-1))+c3*0.5*(x(j)+x(j-1))+a2*(1-exp(-vc*0.5*(z(j)+z(j-1))))-r);
k3=(lambda1(j)-h2*k2)*(c1*0.5*(x(j)+x(j-1))-p*a*0.5*(x(j)+x(j-1))/(a+0.5*(y(j)+y(j-1)))^2)+(lambda2(j)-h2*k2)*(2*r*b*0.5*(y(j)+y(j-1))+c3*0.5*(x(j)+x(j-1))+a2*(1-exp(-vc*0.5*(z(j)+z(j-1))))-r);
k4=(lambda1(j)-h*k3)*(c1*0.5*x(j-1)-p*a*0.5*x(j-1)/(a+0.5*y(j-1))^2)+(lambda2(j)-h*k3)*(2*r*b*0.5*y(j-1)+c3*0.5*x(j-1)+a2*(1-exp(-vc*0.5*z(j-1)))-r);
lambda2(j-1)=lambda2(j)-(h/6)*(k1+2*k2+2*k3+k4);
end
for i=1:N
j=N+2-i;
k1=lambda1(j)*(a1*vc*x(j)*exp(-vc*z(j)))+lambda2(j)*(a2*vc*y(j)*exp(-vc*z(j)))+lambda3(j)*d2;%伴随变量的微分方程，伴随公式
k2=(lambda1(j)-h2*k1)*(a1*vc*0.5*(x(j)+x(j-1))*exp(-vc*0.5*(z(j)+z(j-1))))+(lambda2(j)-h2*k1)*(a2*vc*0.5*(y(j)+y(j-1))*exp(-vc*0.5*(z(j)+z(j-1))))+(lambda3(j)-h2*k1)*d2;
k3=(lambda1(j)-h2*k2)*(a1*vc*0.5*(x(j)+x(j-1))*exp(-vc*0.5*(z(j)+z(j-1))))+(lambda2(j)-h2*k2)*(a2*vc*0.5*(y(j)+y(j-1))*exp(-vc*0.5*(z(j)+z(j-1))))+(lambda3(j)-h2*k2)*d2;
k4=(lambda1(j)-h*k3)*(a1*vc*0.5*x(j-1)*exp(-vc*0.5*z(j-1)))+(lambda2(j)-h*k3)*(a2*vc*0.5*y(j-1)*exp(-vc*0.5*z(j-1)))+(lambda3(j)-h*k3)*d2;
lambda3(j-1)=lambda3(j)-(h/6)*(k1+2*k2+2*k3+k4);
end
u1=max{min{-lambda1/500,0},0};%控制公式
v1=max{min{-lambda3/500,0},0};
u=0.5*(ul+oldu);v=0.5*(vl+oldv);
templ=delta*sum(abs(u))-sum(abs(oldu-u));%对三个变量进行收敛测试
temp2=delta*sum(abs(x))-sum(abs(oldx-x));
temp3=delta*sum(abs(lambda1))-sum(abs(oldlambda1-lambda1));
temp4=delta*sum(abs(lambda2))-sum(abs(oldlambda2-lambda2));
temp5=delta*sum(abs(lambda3))-sum(abs(oldlambda3-lambda3));
temp6=delta*sum(abs(y))-sum(abs(oldy-y));
test=min(min(templ,temp4),min(temp2,temp3),min(temp5,temp6));
end
y(1,:)=x;
y(2,:)=y;
y(3,:)=z;
y(4,:)=u;
y(5,:)=v;
plot (t,y(1,:),'g--');
hold on;
plot (t,y(2,:),'r--');
hold on;
plot(t,y(3, :),'*--');
hold on;
plot(t,y(4,:),'b--');
hold on;
plot(t,y(5,:),'y--');
hold on;
