function root=HalfInterval(f,a,b,eps)
% 二分法求函數(shù)f在[a,b]上的一個零點
% 函數(shù)名：f
% 區(qū)間左端點a
% 區(qū)間右端點b
% 根的精度:eps
% 求出的函數(shù)零點：root
if(nargin==3)
    eps=1.0e-4;
end

f1=subs(sym(f),findsym(sym(f)),a);
f2=subs(sym(f),findsym(sym(f)),b);
if(f1==0)
    root=a;
end
if(f2==0)
    root=b;
end

if(f1*f2>0)
    disp('兩端點函數(shù)值乘積大于0!');
    return;
else
    root=FindRoots(f,a,b,eps);
end

function r=FindRoots(f,a,b,eps)
f_1=subs(sym(f),findsym(sym(f)),a);
f_2=subs(sym(f),findsym(sym(f)),b);
mf=subs(sym(f),findsym(sym(f)),(a+b)/2);
if(f_1*mf>0)
    t=(a+b)/2;
    r=FindRoots(f,t,b,eps);
else if(f_1*mf==0),r=(a+b)/2;
    else if(abs(b-a)<=eps),r=(b+3*a)/4;
        else s=(a+b)/2; r=FindRoots(f,a,s,eps);
        end,end,end