clear all
clear yy
clf

T=10;
xo=.1;
xT=xo/(xo+(1-xo)*exp(-T));

for m=1:10
    Nm=2^m;
    x(1)=xo;
    y(1)=xo;
    dt=T/Nm;
    yy(1,m)=xo;
    for i=2:Nm+1
        a=x(i-1);
        x(i)=a+dt*a*(1-a);
        
        b=y(i-1);
        k1=dt*b*(1-b);
        k2=dt*(b+k1/2)*(1-b-k1/2);
        k3=dt*(b+k2/2)*(1-b-k2/2);
        k4=dt*(b+k3)*(1-b-k3);
        y(i)=b+(k1+2*k2+2*k3+k4)/6;
         yy(i,m)=y(i);
    end
    
    no(m)=Nm;
    akrE(m)=abs((x(Nm+1)-xT)/xT);
    akrRK(m)=abs((y(Nm+1)-xT)/xT);
   
end
loglog(no,akrE,'b',no,akrRK,'r',no,1./no,'--b',no,1./no.^4,'--r')
        for m=1:10
        plot((0:2^m)*T/(2^m),yy(1:2^m+1,m))
        hold on
        end