Метод стрельбы хромающий

[Наталья20]

Всем привет Мне необходимо реализовать метод стрельбы для решения ду 2го порядка, программу я написала на с++(хотя с прогой у меня всегда было туго), дело в том что работает она с ужаснейшей точностью. Ребята, может есть среди вас кто-то, кто сможет подсказать как улучшить точность или же найти возможные ошибки??

Код:
// метод стрельбы.c
#include "stdio.h"
#include "math.h"
#include "iostream"
using namespace std;
const int n = 10000;
// разбиваю ду 2го порядка на систему 1го порядка
double func2(double x, double y, double z)
{
return 6*x;
}
double func1(double x, double y, double z)
{
return z;
}

// Рунге-Кутт для приближения функции

double RK (double y0, double z0, double a, double b, bool met)
{
double k1, k2, k3, k4;
double l1,l2,l3,l4;
double h;
int i;
h = (b-a)/n; // шаг
double x = a; // начальное значение аргумента
double y = y0; // начальное значение функции
double z = z0; // начальное значение производной

for (i = 0; i < n+1; i++)
{

//if (met == 1) cout << x <<" "<< y << " "<< z << endl;

k1 = h*func1(x,y,z);
l1 = h*func2(x,y,z);

k2 = h*func1(x+h/2,y+k1/2,z+l1/2);
l2 = h*func2(x+h/2,y+k1/2,z+l1/2);

k3 = h*func1(x+h/2,y+k2/2,z+l2/2);
l3 = h*func2(x+h/2,y+k2/2,z+l2/2);

k4 = h*func1(x+h,y+k3,z+l3);
l4 = h*func2(x+h,y+k3,z+l3);

x += h;
y += (k1+2*k2+2*k3+k4)/6;
z += (l1+2*l2+2*l3+l4)/6;

}

return y;

}

int main()
{
freopen("data.txt", "r", stdin);//нач приближение считываю с файлика
freopen("result.txt", "w", stdout);

double a = 0.0;// левая граница отрезка
double b = 1.0;// правая граница отрезка

double ya = 0.0;//краевое условие на левой границе
double yb = 1.0;//краевое условие на правой границе

double alpha0;//начальное приближение уравнения(имеется ввиду конечное
//алгебраическое уравнение для решения задачи методом стрельбы)

double alpha1 = 0.0;
double Falpha0, dFalpha0;

double eps = 0.000001, delta, dh = 0.001;
bool met = 0;

cin>>alpha0;
// Метод стрельбы соответственно
do
{

Falpha0 = RK(ya,alpha0,a,b,met)-yb;

dFalpha0 = (RK(ya,(alpha0+eps),a,b,met)-yb-Falpha0)/eps;

alpha1 = alpha0 - Falpha0/dFalpha0;

alpha0 = alpha1;

delta = fabs(Falpha0);
}
while (delta > eps);

met = 1;


cout <<"alpha = "<< alpha0 <<" Runge(alpha0) = "<< RK(ya,alpha0,a,b,met) <<"\n";


return 0;
}