математика — C ++: простой калькулятор четвертого порядка Рунге-Кутта

Я делаю калькулятор для решения нескольких уравнений Рунге-Кутты четвертого порядка, которые мы делали в классе, и хотя я смог заставить калькулятор работать и работать, значения, которые он мне дает, не совсем верны, и я не могу понять, почему.

#include <iostream>
#include <cmath>
#include <iomanip>
#include <cctype>
#include <string>
#include <windows.h>
#include <iostream>
#include <cstdlib>
#include <fstream>
#include <sstream>
#include <vector>
#include <numeric>
#include <algorithm>
using namespace std;
double PI = acos(-1.0);
long double Ks(int,int,long double,long double,long double,long double,long double);

int main ()
{
char answer;
do
int Ans,PoN;
long double X,Y,A,B,h,Fin,K;
cout << "Please enter X and Y used for problem [X Y]: ";
cin >> X >> Y;
cout << "What is the step size? ";
cin >> h;
cout << "What are we approximating to? Y(T) where T = ";
cin >> Fin;
cout << "Which equation form are you using (Please use number)?" << endl;
cout << "1) Ax (+/-) By" << endl;
cout << "2) Ax (+/-) By^2" << endl;
cout << "3) Ay (+/-) By^2" << endl;
cout << "Answer: ";
cin >> Ans;
cout << "You picked: ";
switch(Ans)
{
case 1:
cout << "Ax (+/-) By" << endl;
cout << "Enter A: ";
cin >> A;
cout << "Is the sign in the middle positive(1) or negative(2)? ";
cin >> PoN;
cout << "Enter B: ";
cin >> B;
break;
case 2:
cout << "Ax (+/-) By^2" << endl;
cout << "Enter A: ";
cin >> A;
cout << "Is the sign in the middle positive(1) or negative(2)? ";
cin >> PoN;
cout << "Enter B: ";
cin >> B;
break;
case 3:
cout << "Ay (+/-) By^2" << endl;
cout << "Enter A: ";
cin >> A;
cout << "Is the sign in the middle positive(1) or negative(2)? ";
cin >> PoN;
cout << "Enter B: ";
cin >> B;
break;
}
//cout << fixed << setprecision(8);
cout << " Y[#] |  New Y  | Exact" << endl;
cout << "------|---------|------" << endl;
for(long double i = X+h; i <= Fin; i += h)
{
K = Ks(PoN,Ans,A,B,X,Y,h);
Y = Y + (h/6)*K;
cout << setw(2) << "Y[" << setw(3) << i << setw(1) << "]" << "|" << Y << endl;
}cout << "Would you like to repeat this program?(Y/N) ";
cin >> answer;
}
while (answer == 'y' || answer == 'Y');
if (answer == 'N' || answer == 'n')
{
cout << "Goodbye" << endl;
system("start notepad.exe");
}
system("pause");
return 0;
}
long double Ks(int PoN, int Ans, long double A,long double B,long double X,long double Y,long double h)
{
long double K1=0,K2=0,K3=0,K4=0,K=0;
switch(Ans)
{
case 1:
{
switch(PoN)
{
case 1: //Ax + By
K1 = A*X+B*Y;
K2 = A*(X+((1.0/2)*h)) + B*(Y+(1.0/2)*h*K1);
K3 = A*(X+((1.0/2)*h)) + B*(Y+(1.0/2)*h*K2);
K4 = A*(X+h) + B*(Y+h*K3);
break;
case 2: //Ax - By
K1 = (A*X)-(B*Y);
K2 = A*(X+((1.0/2)*h)) - B*(Y+(1.0/2)*h*K1);
K3 = A*(X+((1.0/2)*h)) - B*(Y+(1.0/2)*h*K2);
K4 = A*(X+h) - B*(Y+h*K3);
//cout << K1 << " " << K2 << " " << K3 << " " << K4 << endl;
break;
}
break;
}
case 2:
{
switch(PoN)
{
case 1: // Ax + By^2
K1 = A*X+pow(B*Y,2);
K2 = A*(X+((1.0/2)*h)) + pow(B*(Y+(1.0/2)*h*K1),2);
K3 = A*(X+((1.0/2)*h)) + pow(B*(Y+(1.0/2)*h*K2),2);
K4 = A*(X+h) + pow(B*(Y+h*K3),2);
break;
case 2: // Ax - By^2
K1 = A*X - pow(B*Y,2);
K2 = A*(X+((1.0/2)*h)) - pow(B*(Y+(1.0/2)*h*K1),2);
K3 = A*(X+((1.0/2)*h)) - pow(B*(Y+(1.0/2)*h*K2),2);
K4 = A*(X+h) - pow(B*(Y+h*K3),2);
break;
}
break;
}

case 3:
{
switch(PoN)
{
case 1: //Ay + By^2
K1 = A*X+B*Y;
K2 = A*(Y+((1.0/2)*h)) + B*(Y+(1.0/2)*h*K1);
K3 = A*(Y+((1.0/2)*h)) + B*(Y+(1.0/2)*h*K2);
K4 = A*(Y+h) + B*(Y+h*K3);
break;
case 2: //Ay - By^2
K1 = A*X-B*Y;
K2 = A*(Y+((1.0/2)*h)) - B*(Y+(1.0/2)*h*K1);
K3 = A*(Y+((1.0/2)*h)) - B*(Y+(1.0/2)*h*K2);
K4 = A*(Y+h) - B*(Y+h*K3);
break;
}
break;
}
}
K = (K1 + 2.0*K2 + 2.0*K3 + K4);
return K;
}

Теперь, чтобы проверить это, я использовал эту проблему:

Где h — размер шага, поскольку он создает новое значение Y для каждой итерации, одна ошибка вызывает ответ на снежный ком после первой итерации. Я подумал, что после десятичной запятой можно рассчитать недостаточно очков, но, как вы видите, я использовал длинные двойные безрезультатно. Использование команды setprecision не работает, она просто добавляет столько нулей на конце, чтобы технически удовлетворить команду.

На этой картинке у меня есть калькулятор Рунге-Кутты, который я нашел слева белым, а моя собственная программа — слева. Как видите, первая итерация верна, но после этого они все больше и больше ошибаются.