Вычислить функцию. Ряды Тейлора

#include <iostream>
#include <cmath>
#include <iomanip>
#define EPS 0.001
using namespace std;
 
double f(double x);
 
int main(){
    
    double x0;
    do{
        cout << "Ââåäèòå x0 = ";
        cin >> x0;
    }while(fabs(x0) <= 1);  
    
    double xk;
    do{
        cout << "Ââåäèòå xk = ";
        cin >> xk;
    }while(fabs(xk) <= 1 || xk < x0 || xk * x0 < 0);
    
    double dx;
    do{
        cout << "Ââåäèòå dx = ";
        cin >> dx;
    }while(dx <= 0);
    
    for (double x = x0; x < xk; x += dx){
        cout << setw(4) << x 
              << setw(12) << setprecision(4) << f(x) 
              << setw(12) << setprecision(4) << log((x+1)/(x-1))
              << setw(12) << setprecision(4) << f(x) - log((x+1)/(x-1)) << endl;
    }
    
    return 0;
}
 
double f(double x){
    double pr = 1 / x;
    double sum = pr;
    int k = 1;
    do{     
        pr *= (2 * k - 1.) / (2 * k + 1.) / x / x;
        sum += pr;
        ++k; 
    }while(fabs(pr) > EPS);
    
    sum *= 2;
    
    return sum; 
}