function ERK_4(a,a0)
%%%%%%%%%%% 比例型延遲微分方程的指數(shù)龍格庫塔方法程序

warning off;
c(1)=0;c(2)=1/2;c(3)=1;
stepnum =50;
q=1/10;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% 變步長格式
m=10;
k = -m:stepnum;
d=1; 
t=zeros(d,m+1+stepnum);
for i=1:m+1
    t(i)=[90*((i-1)/m)+10];
end
for j=1:stepnum
    n= j+m+1;
    t(n)=[(1/q)*t(n-m)];
end

for j = 1:m+stepnum
   h(j) = t(j+1)-t(j);
end
Is=eye(3);
%a=-10;a0=-5;
%%%%%%%%%%%%%%   求拉格朗日插值多項式
syms u;
for i=1:3
    p=1;
    for j=1:3
        n=m+i;
        if j~=i
       p=p*((u/h(n))-c(j))/(c(i)-c(j));
    end
   end
   l(i)=p;
end

for i=1:3
    n=m+i;
    B(i)=(1/h(n))*int(exp((h(n)-u)*a)*l(i),u,0,h(n));
end
for i=1:3
    for j=1:3
        n=m+i;
       A(i,j)=(1/h(n))*int(exp((c(i)*h(n)-u)*a)*l(j),u,0,c(i)*h(n));
    end
end
A=numeric(A);
B=numeric(B);
for i=1:3
    n=m+i;
C(i)=exp(c(i)*(h(n))*a);
end
for n=1:m+1
    k=n-(m+1);
    for i=1:3
    t=k*h(n)+c(i)*h(n);
    Y(n,i)=sin(t);
    y(n)=sin(k*h(n));
end
end
%%%%%%%%%%% 比例型延遲微分方程的指數(shù)龍格庫塔方法
for j=1:stepnum
    n=m+j;
    for i=1:3
  Y(n+1,i)=exp(C(i)*h(n)*a)*y(n)+(1/h(n))*A(i,1)*a0*Y(n-m,1)+(1/h(n))*A(i,2)*a0*Y(n-m,2)+(1/h(n))*A(i,3)*a0*Y(n-m,3);
  y(n+1)=exp(h(n)*a)*y(n)+(1/h(n))*B(1)*a0*Y(n-m,1)+(1/h(n))*B(2)*a0*Y(n-m,2)+(1/h(n))*B(3)*a0*Y(n-m,3);
 
end
 end 
 R=exp(h(n)*a)+h(n)*B*a0*(Is-h(n)*a0*A)^(-1)*C'
t=[-m:stepnum];
plot(t,y,'k');
title('a= -10,a0= 15,初值為y(t)=sin(t)');