Диф уравнения методом прогонки

[АндрейА]

Собственно решение ДУ методом разностных уравнений (метод прогонки)

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

#include <iostream>
#include <math.h>

using namespace std;
const double Cs=1e+20; //начальные условия
const double D=1e-12; //определяют констаты 0,009 и 1,018
const double t=3600;
const double d=0.001;


double Mas[10][5]; //10 строк, 5 столбцов. t/10, d/5
double Alfa[10][5];
double Beta[10][5];

int main()
{
int j;
int i;
int x;
int y;
for (x=0; x<=10; x++) //задает начальные условия
{
y=0;
Mas[x][y]=Cs;
}

//начинаем вычисления коэффициентов прогонки альфа
for (x=0; x<=10; x++)
{
for(y=1; y<=5; y++)
{
Alfa[x][y]=0.009/(1.018-Alfa[x][y-1]*0.009);
}
}
for (x=0; x<=10; x++)
{
for(y=0; y<=5; y++)
{
cout<<" Alfa[i"<<x<<"j"<<y<<"]="<<Alfa[x][y];}
cout <<"\n" ;

}
//начинаем вычисления коэффициентов прогонки бета. Отсюда начинаются проблемы. 4 вложенных цикла не корректно работают друг с другом.
for (i=0; i<=10; i++) //определение начальных условий бета
{
j=0;
Beta[i][j]=Cs;
}
for (x=1; x<=10; x++) //проблемный цикл. На каждом проходе он полностью меняет все значения беты. Должен - закончив расчет одной строки бета и мас - более их не менять.
{
for (y=4; y>=0; y--)
{
for (i=1; i<=10; i++)
{
for(j=1; j<=5; j++)
{
Beta[i][j]=(0.009*Beta[i][j-1]+Mas[x-1][y])/(1.018-Alfa[x][y-1]*0.009);
Mas[x][y]=Alfa[x][y]*Mas[x][y+1]+Beta[i-1][j-1];
}
}
}
}


for (x=0; x<=10; x++) //результаты
{
for(y=0; y<=5; y++)
{
cout<<" "<<Mas[x][y];}
cout <<"\n" ;

}


while(1);
}


Программа работает, но неверно (не может быть отрицательных коэффициентов и превышающих значение Cs).

Вопрос - как можно ещё реализовать совместный расчет второй и третьей матрицы, учитывая что они зависят друг от друга