Rozwiazywanie rownan metoda wielopunktowa

Witam. Mam do zrobienie projekt który wymaga zamiany metody newtona na metode wielopunktowa.
Dokladna tresc zadania:

<url>http://math.uni.lodz.pl/~smietan/pwsz/lab2.pdf<url>. Mi przypadl wzor 3 oraz literka L(na 3 stronie).

Metoda newtona wyglada w ten o to sposob:

#include <iostream>
#include <cmath>

using namespace std;

double f (double x)
{
return(x+0.2)*((x-0.3)*(x-0.3)*(x-0.3));
}

double fp (double x)
{
return (0.3-x)*(0.3-x)* (0.3 +4*x);
}
int main()

{
double x, x0, e, d, k=0;
cin >> x0 >> e;
do {
x= x0 - f(x0)/fp(x0);
k++;
d=fabs(x-x0);
x0=x;
cout <<x << " kroki = " << k << endl;
} while (d>=e);

return 0;
}

f -> funkcja sprawdzajaca metode
x=x0-f(x0)/fp(x0) -> metoda newtona na oszacowanie kolejnego przyblizenia.

Mój problem polega na tym, ze nie jestem do konca pewien jak przerobic to na metode wielopunktowa. Zrobilem wstepny program, ale zapetla mi sie on olbrzymia ilosc razy. Nie wiem czy to ja zrobilem blad w kodzie, czy to tak ma byc. Byłbym niezmiernie wdzieczny za pomoc.

Moj kod:

#include <iostream>
#include <cmath>

using namespace std;



double f(double x)
{
	return pow(x,3)+(3*pow(x,2))+5*x+3;
}
double fp(double x)
{
	return (pow(x,4)/4) +pow(x,3)+((5*pow(x,2))/2)+3*x;
}

double u(double x)
{
	return f(x)/fp(x);
}

double przyblizenie(double x)
{
	return x-u(x)-(f(x-u(x))/fp(x));
}

double wielopunktowa(double x0, double e)
{
	double N,d,x,krok=0;
	cout << " Podaj maksymalna ilosc iteracji" << endl;
	cin >> N;
	do
	{
		x=przyblizenie(x0);
		krok++;
		x0=x;
		cout.precision(16);
		cout <<x << " kroki = " << krok << endl;
	}
	while(krok<=N);
}

double newton(double x0, double e)
{
double x,d,k=0;
do {
x= x0 - f(x0)/fp(x0);
k++;
d=fabs(x-x0);
x0=x;
cout.precision(16);
cout <<x << " kroki = " << k << endl;
} while (d>=e);
}

int main()
{
double x,y;
cin >> x >> y;
cout << "Przyblizenie metoda Newtona" << endl;
newton(x,y);
cout << "Przyblizenie metoda Wielopunktowa" << endl;
wielopunktowa(x,y);

return 0;
}