Решение уравнения теплопроводности по алгоритму Томаса

Я пытаюсь решить дифференциальное уравнение теплопроводности, используя алгоритм Томаса.

Физическая проблема: У нас есть штекер, температура на левой стороне 0температура правой стороны 1,

Для алгоритма Томаса я написал функцию, которая принимает три QVector а также int значение суммы уравнений.

Это мой код:

#include <QCoreApplication>
#include <QVector>
#include <QDebug>
#include <iostream>
using std::cin;

void enterIn(QVector<float> &Array, int Amount_of_elements){
int transit;
for(int i=0;i<Amount_of_elements;i++){
cin>>transit;
Array.push_back(transit);
}
}

QVector<float> shuttle_method(const QVector<float> &below_main_diagonal,
QVector<float> &main_diagonal,
const QVector<float> &beyond_main_diagonal,
const QVector<float> &free_term,
const int N){
QVector <float> c;
QVector <float> d;

for(int i=0;i<N;i++){
main_diagonal[i]*=(-1);
}

QVector<float> x; //result

c.push_back(beyond_main_diagonal[0]/main_diagonal[0]);
d.push_back(-free_term[0]/main_diagonal[0]);

for(int i=1;i<=N-2;i++){
c.push_back(beyond_main_diagonal[i]/(main_diagonal[i]-below_main_diagonal[i]*c[i-1]));
d.push_back(  (below_main_diagonal[i]*d[i-1]  -  free_term[i])  /  (main_diagonal[i]-  below_main_diagonal[i]*c[i-1])  );
}

x.resize(N);
//qDebug()<<x.size()<<endl;

int n=N-1;
x[n]=(below_main_diagonal[n]*d[n-1]-free_term[n])/(main_diagonal[n]-below_main_diagonal[n]*c[n-1]);

for(int i=n-1;i>=0;i--){
x[i]=c[i]*x[i+1]+d[i];
// qDebug()<<x[i]<<endl;
}

return x;
}

int main()
{
QVector <float> alpha;  // below
QVector <float> beta;   // main diagonal * (-1)
QVector <float> gamma;  // beyond
QVector <float> b;      // free term

QVector<float> T;

int cells_x=40;         //amount of equations
alpha.resize(cells_x);
beta.resize(cells_x);
gamma.resize(cells_x);
b.resize(cells_x);
T.resize(cells_x);

float dt=0.2,h=0.1;

alpha[0]=0;
for(int i=1;i<cells_x;i++){
alpha[i]= -dt/(h*h);
}

for(int i=0;i<cells_x;i++){
beta[i] = (2*dt)/(h*h)+1;
}
for(int i=0;i<cells_x-1;i++){
gamma[i]= -dt/(h*h);
}
gamma[cells_x-1]=0;

qDebug()<<"alpha= "<<endl<<alpha.size()<<alpha<<endl<<"beta = "<<endl<<beta.size()<<beta<<endl<<"gamma= "<<gamma.size()<<gamma<<endl;for(int i=0;i<cells_x-1;i++){
T[i]=0;
}
T[cells_x-1]=1;
qDebug()<<endl<<endl<<T<<endl;

//qDebug()<< shuttle_method(alpha,beta,gamma,b,N);

QVector<float> Tn;
Tn.resize(cells_x);
Tn = shuttle_method(alpha,beta,gamma,T,cells_x);
Tn[0]=0;Tn[cells_x-1]=1;
for(int stepTime = 0; stepTime < 50; stepTime++){
Tn = shuttle_method(alpha,beta,gamma,Tn,cells_x);
Tn[0]=0;
Tn[cells_x-1]=1;
qDebug()<<Tn<<endl;
}
return 0;
}

Моя проблема:
когда я компилирую и запускаю его, я получаю это:

Tn  <20 items>  QVector<float>
0   float
0.000425464 float
0.000664658 float
0.000937085 float
0.00125637  float
0.00163846  float
0.00210249  float
0.00267163  float
0.00337436  float
0.00424581  float
0.00532955  float
0.00667976  float
0.00836396  float
0.0104664   float
0.0130921   float
0.0163724   float
0.0204714   float
0.0255939   float
0.0319961   float

Tn  <20 items>  QVector<float>
0   float
-0.000425464    float
0.000643385 float
-0.000926707    float
0.00120951  float
-0.00161561 float
0.00202056  float
-0.00263167 float
0.00324078  float
-0.00418065 float
0.00511726  float
-0.00657621 float
0.00802998  float
-0.0103034  float
0.0125688   float
-0.0161171  float
0.0196527   float
-0.0251945  float
0.0307164   float
1   float

Tn  <20 items>  QVector<float>
0   float
0.000425464 float
0.000664658 float
0.000937085 float
0.00125637  float
0.00163846  float
0.00210249  float
0.00267163  float
0.00337436  float
0.00424581  float
0.00532955  float
0.00667976  float
0.00836396  float
0.0104664   float
0.0130921   float
0.0163724   float
0.0204714   float
0.0255939   float
0.0319961   float

Tn  <20 items>  QVector<float>
0   float
-0.000425464    float
0.000643385 float
-0.000926707    float
0.00120951  float
-0.00161561 float
0.00202056  float
-0.00263167 float
0.00324078  float
-0.00418065 float
0.00511726  float
-0.00657621 float
0.00802998  float
-0.0103034  float
0.0125688   float
-0.0161171  float
0.0196527   float
-0.0251945  float
0.0307164   float
1   float

Снова и снова в цикле.
Я понятия не имею, почему я получаю это.

Может быть, моя ошибка в назначении Tn результат моего Thomas-метода-функции?
или в реализации метода Томаса? или в граничных условиях?