function ys=dbf(f,a,b,a1fa,beta,h,eps)
ff=@(x,y)[y(2),f(y(1),y(2),x)];
xvalue=a:h:b;
n=length(xvalue)
s0=a-0.01;
x0=[a1fa,s0];
flag=0;
y0=rk4(ff,a,x0,h,a,b);
if abs(y0(1,n)-beta)<=eps
flag=1;
y1=y0;
else
s1=s0+1;
x0=[a1fa,s1];
y1=rk4(ff,a,x0,h,a,b);
if abs(y1(1,n)-beta)<=eps
flag=1;
end
end
if flag~=1
while abs(y1(1,n)-beta)>eps
s2=s1-(y1(1,n)-beta)*(s1-s0)/(y1(1,n)-y0(1,n));
x0=[a1fa,s2];
y2=rk4(ff,a,x0,h,a,b);
s0=s1;
s1=s2;
y0=y1;
y1=y2;
end
end
xvalue=a:h:b;
yvalue=y1(1,:);
ys=[xvalue',yvalue'];
function x=rk4(f,t0,x0,h,a,b)
t=a:h:b;
m=length(t);
t(1)=t0;
x(:,1)=x0;
for i=1:m-1
L1=f(t(i),x(:,i));
L2=f(t(i)+h/2,x(:,i)'+(h/2