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
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]))
        yn[3]=yn[2]
        yn[2]=yn[1]
        yn[1]=yp
    end
    ye=1/(1+100*exp(-T))
    er=abs(yn[1]-ye)
    println(h," ",er," ",erold/er)
    erold=er
end