cálculo de raices con regla falsa (matlab)

Hola tengo un código de regla falsa que se hizo en base a uno que se utilizaba para bisección en matlab, pero hay una parte que parece que no se modificó adecuadamente o no comprendo muy bien su función, anexo código:

    
    format long

    syms x
        
    f=input('Ingrese la función f(x):');
    
    fplot(f), grid on
    %saveas(gcf,'Regla Falsa.fig');
    
    xi=input('Ingrese el valor inferior del intervalo (xi):');
    
    xs=input('Ingrese el valor superior del intervalo (xs):');
          
    fi=eval(subs(f,xi));
    fs=eval(subs(f,xs));
    
    if fi==0
        s=xi;
        E=0;
        fprintf('%f es raíz de la función f(x)',xi)
    
    elseif fs==0
        s=xs;
        E=0;
        fprintf('%f es raíz de la función f(x)',xs)
    
    elseif fs*fi<0
        c=0;
        xm=xi-(fi*(xs-xi)/(fs-fi));
        xmn=xm(c+1);
        fm(c+1)=eval(subs(f,xm));
        fe=fm(c+1);
        n(c+1)=c;
        E(c+1)=tol+1;
        error=E(c+1);
        
        while error>tol && fe~=0 && c<niter
            
            if fi*fe<0
                xs=xm;
                fs=eval(subs(f,xs));
            
            else
                xi=xm;
                fi=eval(subs(f,xi));
            end
            
            xa=xm;
            xm=(xi+xs)/2;
            xmn(c+2)=xm;
            fm(c+2)=eval(subs(f,xm));
            fe=fm(c+2);
            E(c+2)=abs(xm-xa);
            %E(c+2)=abs((xm-xa)/xm);
            error=E(c+2);
            c=c+1;
            n(c+1)=c;
        end
        
        if fe==0
           s=xm;
           fprintf('%f es raíz de la función f(x)',xm) 
        
        elseif error<tol
           s=xm;
           fprintf('%f es una aproximación de una raíz de la función f(x) con una tolerancia de %f',xm,tol)
        
        else 
           s=xm;
           fprintf('Fracasó en %d iteraciones',niter) 
        end
        
    else
       fprintf('El intervalo utilizado es inadecuado')         
    end    

    var_names={'n','xmn','fm','Error'};
    val=table(n',xmn',fm',E','VariableNames',var_names);
    
end

En particular la parte del código donde se calcula xm=(xi+xs)/2; es la que me genera dudas, ya que se supone que eso es bisección pero si modifico esa parte por la que corresponde a regla falsa (como está arriba), el código no converge y si borro esa parte, el código no funciona correctamente.