function erroranalysis(trange,y0,N)
%h=zeros(N);e=zeros(N);ediff=zeros(N-1);
a=trange(1);b=trange(2);
for k=1:N
    h(k)=(b-a)/(2^k);
    l=2^k;
    [teuler,yeuler]=eul(@slopefield2,[a,b],y0,h(k));
    [trk2,yrk2]=rk2(@slopefield2,[a,b],y0,h(k));
    [trk4,yrk4]=rk4(@slopefield2,[a,b],y0,h(k));
    %[texact,yexact]=ode45(@slopefield2,teuler,y0,h(k));
    t=a:0.05:b;
    %z=exactsln(t);
    yexact=exactsln(teuler);
    plot(teuler,yeuler,trk2,yrk2)
    %plot(teuler,yeuler,texact,yexact)
    title(['mesh points=',num2str(l),'.'])
    pause
    eEuler(k)=max(abs(yexact-yeuler));
    eRK2(k)=max(abs(yexact-yrk2));
    eRK4(k)=max(abs(yexact-yrk4));
    %format long
    
end
for k=2:N
    ediffEuler(k)=(log(eEuler(k-1))-log(eEuler(k)))/log(2);
    ediffRK2(k)=(log(eRK2(k-1))-log(eRK2(k)))/log(2);
    ediffRK4(k)=(log(eRK4(k-1))-log(eRK4(k)))/log(2);
    
end
subdiv=-log(h/(b-a))/log(2);
plot(subdiv,eEuler,'.',subdiv,eRK2,'.',subdiv,eRK4,'.')
legend('max error')
pause
plot(subdiv,ediffEuler,'.',subdiv,ediffRK2,'.',subdiv,ediffRK4,'.')
legend('log_2 ratio of errors')

 %---------------
    function dydt=slopefield2(t,y)
        %dydt=y; %[0,2],y(0)=1
        %dydt=y.^2*sin(t); %[0,pi], y(0)=0.3
        dydt=y^2;  %[0,1], y(0)=2
        %dydt=sin(t.*y);  %[0,2], y(0)=3
        %dydt=y.^2-t; %[0,1], y(0)=1
        %dydt=y-y.^2; %[0,t], y(0)=0.2
        
    function z=exactsln(t)
         %z=exp(t);
         %z=3./(7+3*cos(t));
         z=2./(1-2*t);