T=10
t0=0
Y(t)=1/(1+100*exp(-t))
f(t,y)=y*(1-y)
df(t,y)=(1-2*y)*f(t,y)
A1=0
A2=1/3
A0=1-A1-A2
B0=23/12+5*A1/12+A2/3
B1=-4/3+2*A1/3+4*A2/3
B2=5/12-A1/12+A2/3
a1=0
a2=2/5
a0=1-a1-a2
bm1=3/8-a1/24
b0=19/24+13*a1/24+a2/3
b1=-5/24+13*a1/24+4*a2/3
b2=1/24-a1/24+a2/3
erold=1
for j=1:12
    global erold
    m=2^j
    h=(T-t0)/m
    yn=Y.([0,-h,-2*h])
    for n=1:m
        tn=t0+h*n
        yp=A0*yn[1]+A1*yn[2]+A2*yn[3]+
            h*(B0*f(tn-h,yn[1])+B1*f(tn-2*h,yn[2])+B2*f(tn-3*h,yn[3]))
        yc=a0*yn[1]+a1*yn[2]+a2*yn[3]+
            h*(bm1*f(tn,yp)+b0*f(tn-h,yn[1])+b1*f(tn-2*h,yn[2])+
                b2*f(tn-3*h,yn[3]))
        yn[3]=yn[2]
        yn[2]=yn[1]
        yn[1]=yc
    end
    ye=1/(1+100*exp(-T))
    er=abs(yn[1]-ye)
    println(h," ",er," ",erold/er)
    erold=er
end