Метод Эйлера для сложного дифура, или что не так с моей математикой?

// Проверено в Dev-C++ 4.9.9.2
#include <cstdlib>
#include <iostream>
#include <math.h>
#include <fstream>
#include <iomanip>
using namespace std;
//------------------------------------------
// ВАЖНО! Чем меньше глобальных переменных -
// тем лучше, но так как это тестовый учебный
// пример, который нужен, чтобы понять принцип,
// то можно закрыть на это глаза.
// ВАЖНО! Крайне опасны однобуквенные глобальные
// переменные.
//------------------------------------------
const int Nv  = 2,//число переменных (уравнений)
          Nt  = 3;//число различных начальных условий
//------------------------------------------
// Параметры модели
//------------------------------------------
double g      = 9.80665, //ускорение силы тяжести, м/(с*с)
       H      = 0.25,    //высота конуса, м
       rho    = 500.0,   //плотность материала конуса, кг/(м*м*м)
       rholiq = 1000.0,  //плотность жидкости, кг/(м*м*м)
       alpha  = 0.5*60.0,//угол полу-раствора конуса, градусов
       f      = 15.0,    //коэффициент сопротивления, Н*с/(м*м*м)
       v0     = 0.0,     //начальная скорость тела, м/с
       t0     = 0.0,     //начальное время моделирования, с
       tmax   = 10.0,    //конечное время моделирования, с
       x0[Nt] = {0.05,
                 0.1,
                 0.15};//три варианта начальных смещений конуса, м
       
double alphar = M_PI*alpha/180.0;//угол полу-раствора конуса, радиан

double h0     = H*pow((rho/rholiq),(1.0/3.0)),//глубина погружения в положении равновесия
       m      = M_PI*rho*H*H*H*tan(alphar)*tan(alphar)/3.0;//масса конуса
       
double c      = M_PI*tan(alphar)*tan(alphar)/m,
       d      = rholiq*g/3.0,
       e      = f/sin(alphar);
//------------------------------------------
// Параметры численного метода
//------------------------------------------
// Желательно, чтобы шаг по времени мог быть
// представлен степенью двойки, а именно
// в форме N/(2^Z). Где N - натуральное, а
// Z - целое число. Например, 1.0/(2.0*2.0*...*2.0).
// 0.125  0.0625
double dt     = 1.0/(2.0*2.0*2.0*2.0*2.0*2.0*2.0*2.0*2.0*2.0*2.0),//шаг по времени, с
       time_c = 0.0,//текущее время
       X[Nv],//переменные
       F[Nv];//производные
int    points = (int)((tmax - t0)/dt) + 1;
double **result;
//------------------------------------------
// Достаточно распространенная форма представления
// систем уравнений при решении численными методами
//------------------------------------------
void SystemODE(double t,double *X,double *F){
  // t - время (в данную систаму ОДУ не входит)
  // F[0] - dx/dt, X[0] - x
  // F[1] - dy/dt, X[1] - y
  // g = g модели
  // c = pi*tg(alpha)*tg(alpha)/m
  // d = rho(liq)*g/3
  // e = f/sin(alpha)
  F[0] = X[1];// dx/dt = y
  F[1] = -g + c*(h0-X[0])*(h0-X[0]) * (d*(h0-X[0]) - e*X[1]);// dy/dt = F(x,y)
}
//------------------------------------------
void EulerStep(double t,double *X,double *F){
  SystemODE(t,X,F);
  for(int i=0;i<Nv;i++)
    X[i] = X[i] + dt * F[i];
}
//------------------------------------------
void Calc(void){
  //Массив для хранения и печати результатов
  result = new double * [points];
  for(int j=0;j<points;j++)
    result[j] = new double [Nt+1];
  //Цикл по вариантам расчета
  for(int i=0;i<Nt;i++){
    //Начальные условия
    time_c = t0;
    X[0]   = x0[i];
    X[1]   = v0;
    //Цикл по расчетным точкам
    for(int j=0;j<points;j++){
      result[j][0]   = time_c;
      result[j][i+1] = X[0];
      EulerStep(time_c,X,F);
      time_c += dt;
    }
  }
  //Печать результата
  fstream fs;
  fs.open("result.txt",std::fstream::out | std::fstream::trunc);
  int digit = 5;
  //fs << std::scientific;
  fs << std::fixed;
  fs << std::setprecision(digit);
  //cout << std::scientific;
  cout << std::fixed;
  cout << std::setprecision(digit);
  for(int j=0;j<points;j++){
    for(int i=0;i<=Nt;i++){
      cout << result[j][i] << "t";
      fs << result[j][i] << "t";
    }
    cout << "n";
    fs << "n";
  }
  fs.close();
  //Удаление массива результатов
  for(int j=0;j<points;j++){
    delete [] result[j];
    result[j] = NULL;
  }
  delete [] result;
  result = NULL;
}
//------------------------------------------
int main(int argc, char *argv[]){
  Calc();
  system("PAUSE");
  return EXIT_SUCCESS;
}
//------------------------------------------