clear
model=4;
ya=0;
yb=1;
a=0;
b=1;
z0=0.04572743438482;
z=z0;
eps=0.00001
x(1)=ya;
x(2)=z0;
[B,u]=rk_func(x,a,b,model,0);
while abs(u(1)-yb)>eps
[B,u]=rk_func(x,a,b,model,0);
f=abs(u(1)-yb);
h=(-u(1)+yb)/B(1,2)
f1=1;
while f1>f 
   [Bn,un]=rk_func(x+h,a,b,model,0);
   f1=abs(un(1)-yb)
   if f1>f
       h=h/2;
   end
end
z=z+h;
x(2)=z;
end
[B,u]=rk_func(x,a,b,model,1);