function p=bisection(FUN,a,b,TOL,N);
%%  function y=bisection(FUN,a,b,tol,N);
%%
%%  The bisection method is to find a solution to f(x)=0 on the interval [a,b]
%%  where f(a) and f(b) have opposite signs
%%
%%  INPUT:
%%  FUN  - the given function to be solved  which is entered as
%%           a string, enclosed in ' ', with the name the m-file of the 
%%           function.
%%  a - the left end point
%%  b - the right end point
%%  TOL - tolerance
%%  N - maximum number of interations 
%%  OUTPUT:
%%  p - approximat solution p or massage of failure
y1=feval(FUN,a);
y2=feval(FUN,b);
if y1*y2>0
   massage='there is no root for this equation'
   return;
end;   
i=1;
while i<=N
    p=a+(b-a)/2;
    if feval(FUN,p)==0
       break;
    elseif (b-a)/2<TOL
       break;
    end;
    i=i+1;
    y=feval(FUN,p);
    if y*y1>0
       a=p;
       y1=y;
    else
        b=p; 
        y2=y;
    end;
end;
if i==N
   masssage='mathod failed after N iteration'
end;             
           
       
