численные методы — Рунге Кутта 4 в C ++ для системы Лоренца

Может кто-нибудь найти ошибку в моем коде? Это работает нормально в специальных точках, но допуск не подходит для метода. Мой Errorstepper довольно прост, поэтому я не думаю, что проблема в этом. Пожалуйста помоги.

#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <conio.h>
#include <time.h>
#include <iostream>
#include <string>
using namespace std;
const int q=10;
const double b=8/3;
double r;
double x,y,z;
struct delta
{
double dx;
double dy;
double dz;
};

double f1(double x , double y)
{
return (q*(y-x));
}
double f2(double x,double y, double z)
{
return ((-x)*z+r*x-y);
}
double f3(double x,double y,double z)
{
return (x*y- b*z);
}
delta calc(double x,double y,double z,double dh=0.1)
{
double k1,k2,k3,k4,m1,m2,m3,m4,p1,p2,p3,p4,dx,dy,dz,a,b,c,h;

h=dh;
k1=h*f1(x,y);
m1=h*f2(x,y,z);
p1=h*f3(x,y,z);

//h+=dh/2;
k2=h*f1((x+k1/2.),(y+m1/2.));
m2=h*f2((x+k1/2.) ,(y+m1/2.),(z+p1/2.));
p2=h*f3((x+k1/2.),(y+m1/2.),(z+p1/2.));

//  h+=dh/2;
k3=h*f1((x+k2/2.),(y+m2/2.));
m3=h*f2((x+k2/2.),(y+m2/2.),(z+p2/2.));
p3=h*f3((x+k2/2.),(y+m2/2.),(z+p2/2.));

//  h+=dh;
k4=h*f1((x+k3),(y+m3));
m4=h*f2((x+k3),(y+m3),(z+p3));
p4=h*f3((x+k3),(y+m3),(z+p3));

dx=(k1+2*k2+2*k3+k4)/6.;
dy=(m1+2*m2+2*m3+m4)/6.;
dz=(p1+2*p2+2*p3+p4)/6.;

a=x+dx;
b=y+dy;
c=z+dz;
delta add;

add.dx=a;
add.dy=b;
add.dz=c;

return add;

}
double Error(double x,double y ,double z, double h)
{
double xi,e,a,b,c;
delta end1=calc(x,y,z,h);
xi=end1.dx;
end1=calc(x,y,z,h/2);
a=end1.dx;
b=end1.dy;
c=end1.dz;
end1=calc(a,b,c,h/2);
e=fabs((xi-end1.dx)/15);
return e;

}
int main ()
{double dx,dy,dz,h,ei,xi,zi,yi,e,t=0,T=0.08,a,b,c,k=0,E,h1,e1,m,E1;
int n=0;
FILE *fff;
fff=fopen("data.txt","w");
cout <<"x0=";
cin>>x;
cout<<"y0=";
cin>>y;
cout<<"z0=";
cin>>z;
cout<<"h=";
cin>>h;
cout<<"r=";
cin>>r;
cout<<"Time=";
cin>>T;
cout<<"Toleranse=";
cin>>E;
xi=x;
double hmin=0.00005;
if (E<=0.0000001)
hmin=hmin/100;
else cout<<"ok"<<endl;

E1=E/10;
do
{e=Error(x,y,z,h);while (e<E1 || e>E)
{
if (e<E1)
h+=hmin;
else if (e>E)
h-=hmin;
else cout<<" ok"<<endl;
e=Error(x,y,z,h);

};/*
while(e>E)
{
h=h/2;
e=Error(x,y,z,h);
};*/
xi=x;
yi=y;
zi=z;
ei=e;

delta konec1=calc(x,y,z,h);

x=end1.dx;
y=end1.dy;
z=end1.dz;

t+=h;
k+=1;
//  cout<<"x="<<x<<" y= "<<y<<" z="<<z<<" Pogreshnost= "<<e<<" Time="<<t<<endl;
fprintf(fff,"%lf, %lf, %lf, %.10lf  ,%lf ,%lf\n", x, y,z,e,t, h);

}
while (t<=T);
fprintf(fff,"Step = %lf",T/k);
fclose(fff);
}