Метод Конечных Разностей для уравнения в частных производных

[illuminates]

Здравствуйте, требуется ваша помощь в нахождении ошибки в программе. Вычисляю по явной схеме. Выдает какие-то очень большие числа.

Задача такая:
[LATEX]U_t=3(1,1-0,5x)U_{xx}+e^t-1[/LATEX]
[LATEX] U(0,t)=0[/LATEX]
[LATEX] U(1,t)=0[/LATEX]
[LATEX]U(x,0)=0.01(1-x)x[/LATEX]

Решение нужно найти с точностью [LATEX]0.0001[/LATEX] на отрезке [LATEX]T=1/a^*, где a^*=\max a(x,t)[/LATEX]
Построить графики функций [LATEX]u(x^*,t), u(x,jt^*)[/LATEX] где [LATEX]x^*=0.6, t^*=T/10, j=1,2,4[/LATEX]

явная разностная схема такая:
([LATEX]$u_t^{j+1}-u_i^j)/\tau=3(1,1-0,5x_i)(u_{i+1}^{j}-2u_i^j+u_{i-1}^j)/h^2+e^{t_j}+1[/LATEX]

код программы:

int main ( void )
{
	setlocale(LC_ALL, "rus");

	int I = 10, J = 30, i, j;
	double  T = 1.0/ pow(3.3, 0.5), h_x = 1.0/ I, h_t = T/ J, epsilon = h_t + pow(h_x, 2), c;
	double **u = new double *[I + 1];
	for (i = 0; i <= I; i++) u[i] = new double [J + 1]; 

	cout<< "Схема может быть неустойчива  при значениях Х :\n";
	for (i = 0; i <= I; i++)
	{
		c = 3 * (1.1 - 0.5 * h_x * i) * h_t * pow(h_x, -2);
		if (c < 0.5) cout << i * h_x << "   ";
	}
	cout <<"\n";
	
	//нулевой слой (j = 0)
 	for (i = 0; i <= I; i++)
	{
		u [i][0] = 0.01 * (1 - i * h_x) * i * h_x;
		//u [i][0] = 1 - i * h_x;  //НУ, несоответствие ГУ и НУ!
	}
		//последующие слои
	for (j = 0; j <= J; j++)
	{
		for (i = 1; i < I; i++) //расчёт j + 1 - го слоя по j-му
		{
			u [0][j + 1] = 0; //ГУ u [0][j + 1] = 1;
			u [I][j + 1] = 0; //ГУ
			u [i][j + 1] = u [i][j] + h_t * (3 * (1.1 - 0.5 * h_x * i) * (u [i + 1][j] -2 * u [i][j] + u [i - 1][j])/ pow(h_x, 2) + exp(h_t * j) - 1);
		}
	}
	int Jv = J/10;

	ofstream out;
    out.open ("D:\\proga7.txt");

	out << "U = U(0.6, t):\n";
	cout << "U = U(0.6, t):\n";
	for (i = 0; i <= J; i++)
	{
		out << h_t * i <<"\t"<< u [6][i] <<"\n";
	    cout << h_t * i <<"\t"<< u [6][i] <<"\n";
	}
	out << "\n U = U(x, 0.33):\n";
	cout << "\n U = U(x, 0.33):\n";
	for (i = 0; i <= I; i++)
	{
		out << h_x * i <<"\t"<< u [i][Jv] <<"\n";
		cout << h_x * i <<"\t"<< u [i][Jv] <<"\n";
	}
	out << "\n U = U(x, 0.66):\n";
	cout << "\n U = U(x, 0.66):\n";
	for (i = 0; i <= I; i++)
	{
		out << h_x * i <<"\t"<< u [i][Jv * 2] <<"\n";
	    cout << h_x * i <<"\t"<< u [i][Jv * 2] <<"\n";
	}
	out << "\n U = U(x, 1.32):\n";
	cout << "\n U = U(x, 1.32):\n";
	for (i = 0; i <= I; i++)
	{
		out << h_x * i <<"\t"<< u [i][Jv * 4] <<"\n";
	    cout << h_x * i <<"\t"<< u [i][Jv * 4] <<"\n";
	}
	out.close();
	getch();
	return 0;
}

выводит следующее: