Решить уравнение при помощи методов Эйлера и Рунге Кутта

[Alexandrietz]

Здравсвуйте, дамы и господа! В институте дали такое задание:
Вот условие задачи
Решить методом Эйлера и Рунге-Кутта уравнение, описывающее смещение x электрона в атоме под действием электрического поля монохроматической световой волны
$$\ddot{x} + \frac{1}{\tau}\dot{x} + {\omega }^{2}x + \gamma {x}^{2} = eE_{0}\cos (pt)$$
[math]\omega[/math] - Круговая частота собственных колебаний электрона в атоме
[math]\tau[/math] - Их время затухания
[math]\gamma[/math] - Коэффициент нелинейности возвращающей силы, действующей на электрон со стороны ядра
[math]e, m[/math] - заряд и масса электрона
[math]E_{0}[/math] - амплитудное значение напряженности электрического поля световой волны
[math]p[/math] - круговая частота волны

Написал программу, касающуюся метода Эйлера, но ее фазовый портрет отвратителен! Помогите найти ошибку

[syntax lang="c"][math]#ifndef FUNCTION
#define FUNCTION

double F1(double E0, double p, double t, double tau, double omega, double gamma, double x, double y);

double F2(double E0, double p, double t, double tau, double omega, double gamma, double x, double y);

#endif[/math]

[math]#include "Function.h"
#include <stdio.h>
#include <math.h>

#define q 1,6E-19
#define m 9,1E-27

double F1(double E0, double p, double t, double tau, double omega, double gamma, double x, double y)
{
	// t - время
  	// E0 - амплитудное значение напряженности электрического поля световой волны
  	// р - круговая частота волны
	// tau - время затухания свободных колебаний
	// omega - круговая частота свободных колебаний электрона в атоме
	// gamma - коэффициент нелинейности возвращающей силы
	return (q * E0 * cos(p * t)/m - y/(tau) - (omega * omega * x) - gamma * x * x); // F1(x, y, t)
}

double F2(double E0, double p, double t, double tau, double omega, double gamma, double x, double y)
{
	// t - время
  	// E0 - амплитудное значение напряженности электрического поля световой волны
  	// р - круговая частота волны
	// tau - время затухания свободных колебаний
	// omega - круговая частота свободных колебаний электрона в атоме
	// gamma - коэффициент нелинейности возвращающей силы
	return y; // F2(x, y, t)
}[/math]

[math]#ifndef EULER
#define EULER
#include <stdio.h>

void Euler(FILE *ef, double E0, double p, double t, double tau, double omega, double gamma, const double a, const double b, const double eps);

#endif[/math]

[math]#include "Function.h"
#include "Euler.h"
#include <math.h>

void Euler(FILE *ef, double E0, double p, double t, double tau, double omega, double gamma, const double a, const double b, const double eps)
{
	double x, x2, y, y2, dt;
	x = a;
	y = b;
	dt = 0.01;
	fprintf(ef, "Initial data x(0) = a = %f; y(0) = b = %f; Epsilon = eps = %f;\n", x, y, eps);
	fprintf(ef, "E0 = %f;", E0);
	fprintf(ef, "\n");
	fprintf(ef, "p = %f;", p);
	fprintf(ef, "\n");
	fprintf(ef, "tau = %f;", tau);
	fprintf(ef, "\n");
	fprintf(ef, "omega = %f;", omega);
	fprintf(ef, "\n");
	fprintf(ef, "gamma = %f;", gamma);
	for(int i = 0;i < 10000; ++i)
		{	
			while (fabs((F1(E0, p, t, tau, omega, gamma, x, y) * dt / 2)) > eps || fabs((F2(E0, p, t, tau, omega, gamma, x, y) * dt / 2)) > eps)	dt/=2.0;
			x2 = x + F1(E0, p, t, tau, omega, gamma, x, y) * dt; // ???????? x ? y ? ???? ?????? ?????? ???? ?? ?? ...
			y2 = y + F2(E0, p, t, tau, omega, gamma, x, y) * dt; // ... ??? ? ? ???? ??????
			x = x2;
			y = y2;
			t += dt;
			fprintf(ef, "%f   %f   %f\n", t , x2, y2);
			dt = 0.01;
		}
}[/math]

[math]#include "Euler.h"
#include <stdio.h>

#define q 1,6E-19
#define m 9,1E -27

int main()
{
	double x, y, E0, p, t, tau, omega, gamma, eps;
	// t - âðåìÿ
  	// E0 - àìïëèòóäíîå çíà÷åíèå íàïðÿæåííîñòè ýëåêòðè÷åñêîãî ïîëÿ ñâåòîâîé âîëíû
  	// ð - êðóãîâàÿ ÷àñòîòà âîëíû
	// tau - âðåìÿ çàòóõàíèÿ ñâîáîäíûõ êîëåáàíèé
	// omega - êðóãîâàÿ ÷àñòîòà ñâîáîäíûõ êîëåáàíèé ýëåêòðîíà â àòîìå
	// gamma - êîýôôèöèåíò íåëèíåéíîñòè âîçâðàùàþùåé ñèëû
	// epsilon = eps - 
	printf("Enter epsilon\n");
	scanf("%lf", &eps);
	printf("Enter x, y\n");
	scanf("%lf %lf", &x, &y);
	printf("Enter E0:\n");
	scanf("%lf", &E0);
	printf("Enter p:\n");
	scanf("%lf", &p);
	printf("Enter t:\n");
	scanf("%lf", &t);
	printf("Enter tau:\n");
	scanf("%lf", &tau);
	printf("Enter omega:\n");
	scanf("%lf", &omega);
	printf("Enter gamma:\n");
	scanf("%lf", &gamma);
	printf("\n");
	FILE *ef;
	ef = fopen("Euler.txt", "w");
	Euler(ef, E0, p, t, tau, omega, gamma, x, y, eps);
	fclose(ef);
	return 0;
}[/math][/syntax]