Метод Рунге-Кутта

[Zhiend]

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

#include <stdio.h>
#define _USE_MATH_DEFINES
#include <math.h>
#define N 43200.0

void Compute(double &x0, double y, double &z0, double k1, double k2, double k3, double k4, double h, double i, FILE *stream) {
	double KX1, KX2, KX3, KZ1, KZ2, KZ3;
	double x, z, t;
	x = x0;
	z = z0;
	t = i;
	//Уравнения
	KX1 = -(k1 * x * y) - (k2 * x * z) + (2 * k3 * y) + (k4 * z);
	KZ1 = (k1 * x * y) - (k2 * x * z) - (k4 * z);
	x = x0 + (h / 3) * KX1;
	z = z0 + (h / 3) * KZ1;
	t = i + (h / 3);
	KX2 = -(k1 * x * y) - (k2 * x * z) + (2 * k3 * y) + (k4 * z);
	KZ2 = (k1 * x * y) - (k2 * x * z) - (k4 * z);
	x = x0 + (2 * h / 3) * KX2;
	z = z0 + (2 * h / 3) * KZ2;
	t = i + (2 * h / 3);
	KX3 = -(k1 * x * y) - (k2 * x * z) + (2 * k3 * y) + (k4 * z);
	KZ3 = (k1 * x * y) - (k2 * x * z) - (k4 * z);
	x0 = x0 + (h / 4) * (KX1 + 3 * KX3);
	z0 = z0 + (h / 4) * (KZ1 + 3 * KZ3);
	fprintf(stream, "X = %lf\nZ = %lf\n", x0, z0);
}

int main() {
	FILE* stream;
	fopen_s(&stream, "result", "w");
	//начальные условия задачи Коши
	double x0 = pow(10, 6), y0 = 3.7 * pow(10,16), z0 = pow(10, 12), t = 0.0;
	//Коэффициенты в уравнениях
	double k1 = 1.63 * pow(10, -16), k2 = 4.66 * pow(10, -16);
	// k3 и k4 изменяются дважды в сутки
	double k3_1 = 0, k4_1 = 0, k3_2 = expl(-7.601), k4_2 = expl(-22.62);
	// шаг интегрирования
	double h = 5;
	for (t; t < N; t = t + h) {
		if (t != N / 2) {
			Compute(x0, y0, z0, k1, k2, k3_1, k4_1, h, t, stream);
		}
		else {
			Compute(x0, y0, z0, k1, k2, k3_2, k4_2, h, t, stream);
		}
	}
}