połączenie bisekcji i metody Newtona-Rhapsona

Witam,
mam do napisania program który znajdzie miejsca zerowe funkcji używając do tego bisekcji i metody N-R
na razie mój kod wygląda następująco:

#include <iostream>
#include <cmath>
using namespace std;

double funkcja(double x)
{      
       return pow(x,2) - 2*x - 2;     
}

double p1(double x)
{
       return 2*x - 2;       
}

int main()
{
double xr=0.7, x=0.4, xl=0.4, tol = 0.000001, xm, xi;
int bi=0, nr=0, i=0;

cout.precision(10);

while(abs((xr-xl)/xr) >= tol)
{
if( xl <= ( x - (funkcja(x) / p1(x)) ) && ( x - (funkcja(x) / p1(x) )) <= xr )
//if( ((x-xl)*p1(x)-funkcja(x))>=0 && ((x-xr)*p1(x)-funkcja(x))<=0 )
{
    x = xl;  
    xi = (x - (funkcja(x)/p1(x)));
    for(i = 1; abs( (xi-x) / xi ) >= tol; i++)
    {     
          xl = xi;
          xi = (x - (funkcja(x)/p1(x)));    
    }
    nr++;
}
else
{
    x = xl;
    xm = (xl+xr)/2;
    cout<<xl<<" "<<xr<<" "<<xm<<endl;
    if((funkcja(xl) * funkcja(xm)) < 0)
    {
             xr = xm;
    }
    else 
    {
             xl = xm;  
    }
    
    bi++;
} 
}  
    
cout << x << endl;
cout << "bi = " << bi << endl;
cout << "nr = " << nr << endl;
system("pause");
return 0;    
}

pytanie: gdzie popełniłem błąd ? czy while() powinien wyglądać właśnie tak ? może pomieszałem gdzieś zmienne ?
zwracany wynik nie zgadza się z oczekiwanym
jeżeli ktoś zauważył jakiś błąd, byłbym wdzięczny za jego wskazanie

nie ma miejsca zerowego w przedziale [0.4, 0.7]