[C++] Решение системы диференциальных уровнений методом Рунге-Кутта

[Dp12]

Доброво времени суток. Есть задача составить программу решения системы диференциальных уравнений методом Рунге-Кутта. Когда запуская программу выводить очень больший числа. Не знаю в чем проблема.
Приклепляю код программы.

#include <iostream>
#include <cmath>

using namespace std;

double f1(double x1, double x2)
{
	return 12 - x1 - 36 * x2;
}

double f2(double x3)
{
	return x3;
}

void enterData(double *a, double *b, double *h) {
	cout << "Введите границы отрезка a и b: \n";
	cin >> *a >> *b;
	cout << "\nВведите шаг h:\n";
	cin >> *h;

}

int main(void)
{
	setlocale(LC_ALL, "rus");
	double a, b, h;
	int n;
	double *x, *x1, *x2, *k11, *k12, *k13, *k14, *k21, *k22, *k23, *k24;

	enterData(&a, &b, &h);
	

    a = 0, b = 0.5, h = 0.1;
	//n = (b - a) / h;
	n = (a + b) / h;

	
	x = new double[n + 1];
	x1 = new double[n + 1];
	x2 = new double[n + 1];
	k11 = new double[n + 1];
	k12 = new double[n + 1];
	k13 = new double[n + 1];
	k14 = new double[n + 1];
	k21 = new double[n + 1];
	k22 = new double[n + 1];
	k23 = new double[n + 1];
	k24 = new double[n + 1];

	cout << "\nВведите начальные условия.\nx = ";
	cin >> x1[0];
	//cout << "\ny = ";
	cin >> x2[0];
	cout << endl;
x[0] = 0;
	
	for (int i = 1;i <= n;i++) {
		x[i] = a + i*h;
k11[i] = h*f2(x2[i]);
		k21[i] = h*f1(x1[i], x2[i]);

		k12[i] = h*f2(x2[i] + k21[i] / 2);
		k22[i] = h*f1(x1[i] + k11[i] / 2, x2[i] + k21[i] / 2);

		k13[i] = h*f2(x2[i] + k22[i] / 2);
		k23[i] = h*f1(x1[i] + k12[i] / 2, x2[i] + k22[i] / 2);

		k14[i] = h*f2(x2[i] + k23[i]);
		k24[i] = h*f1(x1[i] + k13[i], x2[i] + k23[i]);

		x1[i + 1] = x1[i] + (k11[i] + 2 * k12[i] + 2 * k13[i] + k14[i]) / 6;
		x2[i + 1] = x2[i] + (k21[i] + 2 * k22[i] + 2 * k23[i] + k24[i]) / 6;

	}
for (int i = 0; i < n + 1; i++) {
		cout << "X[" << i << "] = " << x[i] << " \n";
	}
	for (int i = 0; i < n + 1; i++) {
		cout << "X1[" << i << "] = " << x1[i] << "    X2[" << i << "] = " << x2[i] <<
			";\n";
	}


	return 0;

}

[WorldMaster]

После инициализации массивов обнулите их как положено.

memset(x, 0, sizeof(double)*(n + 1));
	memset(x1, 0, sizeof(double)*(n + 1));
	memset(x2, 0, sizeof(double)*(n + 1));
	memset(k11, 0, sizeof(double)*(n + 1));
	memset(k12, 0, sizeof(double)*(n + 1));
	memset(k13, 0, sizeof(double)*(n + 1));
	memset(k14, 0, sizeof(double)*(n + 1));
	memset(k21, 0, sizeof(double)*(n + 1));
	memset(k22, 0, sizeof(double)*(n + 1));
	memset(k23, 0, sizeof(double)*(n + 1));
	memset(k24, 0, sizeof(double)*(n + 1));