Численное решение задачи Коши

class gridfun
{
private:
    double *x;
    double *y;
    int n;
public:
    gridfun(double, double, double);
    double get_x(int i);
    double get_y(int i);
    int get_n();
    void set_x(int i, double value);
    void set_y(int i, double value);
};

int gridfun::get_n()
{
    return n;
}
 
double gridfun::get_x(int i)
{
    return *(x + i);
}
 
double gridfun::get_y(int i)
{
    return *(y + i);
}
 
void gridfun::set_x(int i, double value)
{
    *(x + i) = value;
}
 
void gridfun::set_y(int i, double value)
{
    *(y + i) = value;
}
 
 
gridfun::gridfun(double h, double a, double b)
{
    n = (int)((b - a) / h);
    if (b - a - n*h > 0) n++;
    x = new double[n];
    y = new double[n];
    for (int i = 0; i < n + 1; i++)
    {
        *(x + i) = a + i*h; *(y + i) = 0;           
        if (i == n) *(x + i) = b;                   
    }
}

// Work_1
// Модуль реализует численное решение уравнения y'=x^2(y^2+1)
double function(double x, double y)             
{
    return x*x*(y*y + 1);
}
 
 
// Solve Work_1
double ExSol(double x)// точное решение уравнения sin((x^3)/3) / cos((x^3)/3)
{
    return sin(x*x*x / 3) / cos(x*x*x / 3);
}

/ Метод Эйлера
void Euler(gridfun &f)
{
    int n = f.get_n();
    double b, c, xt, xp, yp, yt;
    f.set_y(0, .0);
    for (int i = 1; i < n + 1; i++)
    {
        xt = f.get_x(i); xp = f.get_x(i - 1);
        b = -1.0 / ((xt - xp)*xt*xt); yp = f.get_y(i - 1);
        c = 1 + yp / ((xt - xp)*xt*xt);
        yt = 0.5*(-pow(b*b - 4.0*c, 0.5) - b);
        f.set_y(i, yt);
    }
}
 
// Метод Рунге - Кутта 3 порядка
void Runge_Kutta3(gridfun &f)
{
    int n = f.get_n();
    double xt, xp, yp, yt;
    double k1, k2, k3;
    f.set_y(0, .0);
    for (int i = 1; i < n + 1; i++)
    {
        xt = f.get_x(i); xp = f.get_x(i - 1);
        yp = f.get_y(i - 1);
        k1 = (xt - xp)*function(xp, yp);
        k2 = (xt - xp)*function(xp + (xt - xp) / 2, yp + k1 / 2);
        k3 = (xt - xp)*function(xp + (xt - xp), yp + 2 * k2 - k1);
        yt = yp + (k1 + 4 * k2 + k3) / 6;
        f.set_y(i, yt);
    }
}

gridfun *f1, *f2, *f3, *f4, 
*ff1, *ff2, *ff3;
 
int main(void)
{
f1 = new gridfun(h, 0.0, 1.0);
f2 = new gridfun(h, 0.0, 1.0);
 
Euler(*f1);
Runge_Kutta3(*f2);
 
cout << "Таблица значений (неявный метод Эйлера):" << endl;
        cout << " N | x[i]       | y[i]       | y(x[i]) | eps[i]" << endl;
        for (int i = 0; i < f1->get_n() + 1; i++)
        {
            cout << " " << i << " ";
            cout << fixed << setprecision(8) << "| " << f1->get_x(i) << " ";
            cout << "| " << f1->get_y(i) << " ";
            cout << "| " << ExSol(f1->get_x(i)) << " ";
            cout << "| " << fabs(f1->get_y(i) - ExSol(f1->get_x(i))) << endl;
            if (i == f1->get_n())
            {
                cout << "|eps| = " << fabs(f1->get_y(i) - ExSol(f1->get_x(i))) << endl;
            };
        }
}
 
    cout << "\nТаблица значений (метод Рунге-Кутты 3):" << endl; // Метод Рунге - Кутты 3 порядка
        cout << " N | x[i]       | y[i]       | y(x[i]) | eps[i]" << endl;
        for (int i = 0; i < f3->get_n() + 1; i++)
        {
            cout << " " << i << " ";
            cout << fixed << setprecision(8) << "| " << f3->get_x(i) << " ";
            cout << "| " << f3->get_y(i) << " ";
            cout << "| " << ExSol(f3->get_x(i)) << " ";
            cout << "| " << fabs(f3->get_y(i) - ExSol(f3->get_x(i))) << endl;
            if (i == f3->get_n())
            {
                cout << "|eps| = " << fabs(f3->get_y(i) - ExSol(f3->get_x(i))) << endl;
            };
        }

        *(x + i) = a + i*h; *(y + i) = 0;           
        if (i == n) *(x + i) = b;