Метод Ньютона-Рафсона в С++

Для решения алгебраических уравнений используем метод Ньютона (его еще называют итерационным методом Ньютона-Рафсона). Предположим, необходимо найти корень уравнения $f\left(x \right)=0$ в окрестности точки $x=x_{0}$. Последовательный поиск приближений для корня уравнения выполняется по следующей схеме: в точке текущего приближения для корня уравнения проводится касательная к графику функции $f\left(x \right)$, точка пересечения этой касательной с осью абсцисс является новым приближением для корня. Такая последовательность действий сводится к рекуррентной формуле: $$x_{n+1}=x_{n}-f\left(x_{n} \right)/f'\left(x_{n} \right)$$ Здесь расмотрим случай, когда исходная функция представляет собой полином. На входе задаются только коэффициенты этого полинома, набор которых можно представить в виде элементов массива. $$f\left(x \right)=a_{0}+a_{1}x+a_{2}x^{2}+a_{3}x^{3}+a_{4}x^{4}+...$$ Кроме того пользователем вводится и начальное приближение корня. Это означат, что надо иметь представление о том, где приблизительно находится корень. Другими словами его надо отделить некоторым интервалом.

#include 
#include 
using namespace std;
int main(){
//Порядок полинома:
const int N=3;
//Индексные переменные и число итераций:
int i,k,n=10;
//Аргумент, функция и производная:
double x,f,df;
//Массив коэффициентов полинома-функции:
double a[N+1];
//Массив коэффициентов полинома-производной:
double b[N];
//Ввод коэффициентов функции-полинома:
cout << "Function: ";
cin >> a[0];
for(i=1; i < N+1 ;i++){
 cin >> a[i];
 //Вычисление коэффициентов для производной:
 b[i-1]=i*a[i];
}
//Начальное приближение:
cout << "Enter x0 = ";
cin >> x;
//Последовательные итерации:
for(k=1; k <= n; k++){
 f=a[0];
 df=0;
 for(i=1; i < N+1; i++){
 f+=a[i]*pow(x,i);
 df+=b[i-1]*pow(x,i-1);
 }
 x-=f/df;
}
//Результат:
cout << "x = " << x << endl;
cout << "f(" << x << ") = " << f << endl;
return 0;
}