N=20;
f=@(x) (x+0.01*sin(x)^2)^(-1/2);
% f=@(x) x^(-1/2)

for m=1:5

    N=N*2;
    dx=1/N;

    y=zeros(N+1,1);
    x=zeros(N+1,1);
    y(1)=0;
    x(1)=0;

    for n=1:N
        x(n+1)=n*dx;
        y(n+1)=y(n)+dx*f(x(n+1));
    end


    err(m)=abs(y(end)-2);
    Dx(m)=dx;

end
%%
loglog(Dx,err,Dx,Dx.^0.5)
legend('error','half order')
xlabel('dx')
ylabel('error')
