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);