%對于指定的a,畫出真解和近似解,紅色的為真解
%解方程y'=D^(1-a)(-y+g(t))  y(0)=0  其中g(t)=t^2+2t^(2-a)/gamma(3-a);
%真解是y=t^2;
clear
global a
a=0.6;   
n_step=100;
h=1/n_step;
y(1)=0;
for k=1:n_step
    t=(k-1)*h;
    fy(k)=-y(k)+t^2+2*t^(2-a)/gamma(3-a); 
    out=0;
    for j=0:k-1
          out=out+b(j,k-1)*fy(j+1);
    end
    y(k+1)=h^(a-1)/gamma(1+a)*out*h+y(k);
end
x=0:h:1;
plot(x,y);
hold on;
y1=x.^2;
plot(x,y1,'r');
xlabel('t')
ylabel('y')
title('Eular')
text(0.9,y(91),'\leftarrow y_h',...
     'HorizontalAlignment','left')  
text(0.7,0.7^2,'t^2 \rightarrow',...
     'HorizontalAlignment','right')