p=0:-1:-15;
h=2.^p;
f = @(t,y) 5*y-3;
ytrue = @(t) 2/5*exp(5*t)+3/5;
y0=1;
ye=zeros(size(h,2),1/h(1,end)+1);
yr=zeros(size(h,2),1/h(1,end)+1);
ye(:,1)=y0*ones(size(h,2),1);
yr(:,1)=y0*ones(size(h,2),1);

for j=1:size(h,2)
    for i=1:1/h(1,j)
        ye(j,i+1)=ye(j,i)+h(1,j)*f(i*h(1,j),ye(j,i));
        yr(j,i+1)=yr(j,i)+h(1,j)*f(i*h(1,j)+h(1,j)/2,yr(j,i)+h(1,j)/2*f(i*h(1,j),yr(j,i)));
    end
end

for j=1:size(h,2)
    r(1,j)=ye(j,1/h(1,j)+1);
    r2(1,j)=yr(j,1/h(1,j)+1);
end
yref=ytrue(1)*ones(1,size(h,2));
plot(h,abs(r-yref))
hold on
plot(h,abs(r2-yref),'r')
%plot(h,h.^2,'r:')
xlabel('h')

figure
plot(r)
hold on
plot(r2,'r')
plot(yref,'g')
xlabel('i')


%plot(p,ye(:,end),'b')
%hold on
%plot(p,yr(:,end),'r')