//1st fractional step for diffusion term in eta-space
//Implicit method using Thomas algorithm for tri-diagonal matrix

	scalarField a(etamax+1, 0), b(etamax+1, 0), c(etamax+1, 0);

        for(label j=0 ; j <= etamax ; j++)
        {
		if(j == int(Neta[1]) || j == int(Neta[2]) || j == int(Neta[3]))
                {
                	a[j] = -SDRCMC[j+k*(etamax+1)] * runTime.deltaT().value() /(  (etaValue[j] - etaValue[j-1]) * (etaValue[j+1] - etaValue[j-1]) / 2 ) ;
                        b[j] = 1+2*SDRCMC[j+k*(etamax+1)] * runTime.deltaT().value() /(  (etaValue[j+1] - etaValue[j]) * (etaValue[j] - etaValue[j-1]) ) ;
                        c[j] = -SDRCMC[j+k*(etamax+1)] * runTime.deltaT().value() /(  (etaValue[j+1] - etaValue[j]) * (etaValue[j+1] - etaValue[j-1]) / 2 );
			
                }
                else
                {
                	scalar Gam = SDRCMC[j+k*(etamax+1)] * runTime.deltaT().value() / (deta[j] * deta[j]) ;
                        a[j] = -Gam;
                        b[j] = (1+2*Gam);
                        c[j] = -Gam;
                }
        }

        for(label i = 0 ; i<Ysize ; i++)
        {
                scalarField QiCMC_temp(etamax+1,0), QiCMC_temp2(etamax+1, 0);

                for(label j=0 ; j<=etamax ; j++)
                {
                        QiCMC_temp[j] = QiCMC[(k*(etamax+1)+j)*Ysize + i] ;
                }

                        TDMA( a,  b,  c, QiCMC_temp, QiCMC_temp2, etamax+1);

                for(label j=0 ; j <= etamax ; j++)
                {
                        if( j != 0 || j != etamax)
                        {
                         QiCMCNew[(k*(etamax+1)+j)*Ysize + i]  = QiCMC_temp2[j];
                        }
                        else
                        {
                         QiCMCNew[(k*(etamax+1)+j)*Ysize + i]  = QiCMC_temp[j];
                        }
                }
        }