clear
model=2;
ya=1;
yb=2;
a=0;
b=1;
z0=0;
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);
h=(-u(1)+yb)/B(1,2)
z=z+h;
x(2)=z;
end
[B,u]=rk_func(x,a,b,model,1);