分数阶微分方程求解--预估校正法

function [t,y] = fo_solution(h,tfinal,ini_val,alpha,delay)
%%%%%%%%%%%
%%%%%%%%%%%
%	fo_solution.m: 分数阶微分方程求解函数，这个一般不需要改，
%		参数：
%			h(步长),
%			tfinal(终止时间),
%			ini_val(系统初始值),
%%			alpha(阶数),
%			delay(时延)，
%		要求解带有时延的系统，时延只能是离散时延。
%		若求解不带有时延的系统，时延默认为0。
%		返回值：
%			t是时间向量，
%			y是系统的值，第i行是变量i的数值解。
%
%	自己创建equation_func.m: 写系统的微分方程。
%		参数： 
%			t(当前时刻),
%			x(当前时刻系统的值), x为列向量
%			xd(时延项的系统值), xd为列向量
%		返回值：为列向量。
    if nargin<5
        delay=0;
    end
    len = length(ini_val);
    N = tfinal/h;
    t = linspace(0,N,N+1)*h;
    k = delay/h;
    y0 = ini_val;
    y = zeros(len,N);
    f0 = equation_func(0,y0,y0);
    x1=zeros(1,N);
    x2=zeros(1,N);
    x3=zeros(1,N);
    for n = 0:N-1
        disp(['运行次数： ',num2str(n),' t = ',num2str(n*h)]);
        b=@(j)  h^alpha/alpha*(    (n+1-j)^alpha-(n-j)^alpha  );
        s=zeros(len,1);
        for j = 0:n
            if j-k<=0
                x_j_k=y0;
            else
                x_j_k=y(:,j-k);
            end
         if j==0
             x_j=y0;
         else
             x_j=y(:,j);
         end
            s = s+b(j)*equation_func(j*h,x_j,x_j_k);
        end
        
        yp = y0+s/gamma(alpha);
        a=@(j) (n-j+2)^(alpha+1)+(n-j)^(alpha+1)-2*(n-j+1)^(alpha+1);
        S0=(n^(alpha+1)-(n-alpha)*(n+1)^alpha)*f0;
        S1 = zeros(len,1);
        for j = 1:n
            if j-k<=0
                x_j_k=y0;
            else
                x_j_k=y(:,j-k);
            end
            x_j=y(:,j);
            S1 = S1+a(j)*equation_func(j*h,x_j,x_j_k);
        end
        S=S0+S1;
        if n+1-k<=0
            x_n_k=y0;
        else
            x_n_k=y(:,n+1-k);
        end
        y(:,n+1) = y0+h^alpha/gamma(alpha+2)*equation_func(n*h,yp,x_n_k)+S*h^alpha/gamma(alpha+2);
    end
    y = [y0,y];
   
end

有错误之处，请留言。带有时变时延的分数阶微分方程有没有大佬会的，跪求!!! !!!