Добавить вычисление с точносью для метода Рунге-Кутты 2-го порядка

#include <iostream>
#include <cmath>
using namespace std;
float fun(float x, float y)
{
    return x/(2*y);
}
 
void run_kut(float a, float b, float y0, int n)
{
    float x, y, x1, y1, h, k1, k2, k3, k4;
    x = a;
    y = y0;
    h = (b - a) / n;
    for (int i = 0; i<n; i++)
    {
        k1 = h*fun(x, y);
        x1 = x + h / 2.;
        y1 = y + k1 / 2.;
        k2 = h*fun(x1, y1);
        y1 = y + k2 / 2.;
        k3 = h*fun(x1, y1);
        x += h;
        y1 = y + k3;
        k4 = h*fun(x, y1);
        y += (k1 + 2. * k2 + 2. * k3 + k4) / 6.;
        cout << "x[" << i + 1 << "]=" << x << "\t" << "y[" << i + 1 << "]=" << y << "\n";
    }
}
int main()
{
    setlocale(LC_ALL, "Russian");
    int n;
    float a, b, y0,e;
    system("Cls");
    cout << "Ведите a="; cin >> a;
    cout << "Ведите b="; cin >> b;
    cout << "Ведите  y0="; cin >> y0;
    n = (b-a)/pow(0.0001,1./4);
    cout<<"Шаг = " << (b - a) / n<<endl;
    cout << "\nМетод Рунге-Кутты:\n";
    run_kut(a, b, y0, n);
    cout << "\n";
    return 0;
    system("pause");
}

#include <iostream>
using namespace std;
 
double fdY(double X, double Y)
{
    return X / (2 * Y);
}
 
void run_kut(double Y0,double a, double b, int n)
{
    double h = (b - a) / n;
 
    double X = a;
    double Y = Y0;
 
    double k1Y, k2Y, k3Y, k4Y;
 
    cout << "X\tY" << endl;
    do
    {
        cout << X << "\t" << Y << endl;
        k1Y = h * fdY(X, Y);
        k2Y = h * fdY(X + h / 2, Y + k1Y / 2);
        k3Y = h * fdY(X + h / 2, Y + k2Y / 2);
        k4Y = h * fdY(X + h, Y + k3Y);
 
        X += h;
        Y += (k1Y + 2 * k2Y + 2 * k3Y + k4Y) / 6;
    } while (X <= b + h / 2);
}
 
int main()
{
    double Y0 = 1.0;
    double a = 0.0;
    double b = 1.2;
    int n = 10;
    run_kut(Y0, a, b, n);
    system("pause");
    return 0;
}