%%---------------------------------------------------------------------------
function [x,fx,xflag]=find_roots1(f,c1,c2,rel_err)
% f is the name of the function to be solved
% c1,c2 are the min and max bounds of x respectively
% rel_err is the relative error,that isabs(f1n-f2n)/abs(f1-f2)
% where f1n, f2n are the function values at the c1n and c2n
% f1n, f2n are the function values at c1 and c2
% xflag is -1 if there's no root in range (c1,c2)
% xflag equals to 1 if there is a root in range (c1,c2)
% x is the root solved
% fx is its value
% example
% [x,fx,xflag]=find_roots1(@sin,0.1,pi+0.1,1e-6)
%{
c1=0.01;
c2=pi+0.1;

rel_err=1e-3;
f=@sin;
%}
xflag=-1;
f1=feval_r(f,c1);
f2=feval_r(f,c2);
abs_f1_2=abs(f1-f2);
if abs(f1) <= rel_err*abs_f1_2
x=c1;
fx=f1;
xflag=1;
return;
elseif abs(f2) <= rel_err*abs_f1_2
x=c2;
fx=f2;
xflag=1;
return;
elseif f1*f2 >————0
x=nan;
fx=nan;
return;
else
fori1=1:100
c0=(c1+c2)/2;
f0=feval_r(f,c0);
if abs(f0) <= rel_err*abs_f1_2
x=c0;
fx=f0;
xflag=1;
return;
elseif f1*f0 >0
c1=c0;
f1=f0;
else
c2=c0;
f2=f0;
end
if abs(f1-f2)<rel_err*abs_f1_2
xflag=1;
x=c0;
fx=f0;
return;
elseif abs(f1-f2)>100*abs_f1_2
xflag=-1;
x=nan;
fx=nan;
return;
end
end
end
%%---------------------------------------------------------------------------