用指数龙格库塔方法求解时滞（延迟）微分方程！

				function ERK4(h,tao,N)
				%%%%%%%%%%%延迟微分方程的指数龙格库塔方法程序
				
				c(1)=0;c(2)=1/2;c(3)=1;
				m=ceil(tao/h);
				Is=eye(3);
				a=-10;
				a0=-5;
				
				syms u;
				for i=1:3
				    p=1;
				    for j=1:3
				        if j~=i
				       p=p*((u/h)-c(j))/(c(i)-c(j));
				    end
				   end
				   l(i)=p;
				end
				
				q=exp((h-u)*a);
				for i=1:3
				    B(i)=(1/h)*int(q*l(i),u,0,h);
				end
				for i=1:3
				    for j=1:3
				       A(i,j)=(1/h)*int(exp((c(i)*h-u)*a)*l(j),u,0,c(i)*h);
				    end
				end
				A=numeric(A);
				B=numeric(B);
				
				for i=1:3
				C(i)=exp(c(i)*h*a);
				end
				for n=1:m+1
				    k=n-(m+1);
				    for i=1:3
				    t=k*h+c(i)*h;
				    Y(n,i)=sin(t);
				end
				end
				%%%%%%%%%%%延迟微分方程的指数龙格库塔方法
				y(m+1)=0;
				for n=m+1:N
				    for i=1:3
				  Y(n+1,i)=exp(h*a)*y(n)+h*A(i,1)*a0*Y(n-m,1)+h*A(i,2)*a0*Y(n-m,2)+h*A(i,3)*a0*Y(n-m,3);
				  y(n+1)=exp(h*a)*y(n)+h*B(1)*a0*Y(n-m,1)+h*B(2)*a0*Y(n-m,2)+h*B(3)*a0*Y(n-m,3);
				end
				 end 
				 R=exp(h*a)+h*B*a0*(Is-h*a0*A)^(-1)*C'
				t=[1:N+1]*h;
				plot(t,y,'k');
				title('a= -10,a0= -5,初值为y(t)= sin(t)')