#include <iostream>
#include <cmath>
#include <cstdlib>

using namespace std;

void trdiag(int kdim,int n,double *bef,double *diag,double *aft,double *rhsol)
{
		
	double eps=1.e-8;
	for(int i=1;i<n;i++)	
	{
		if(abs(diag[i-1])<eps)
		{
			cout<<"Divide by zero"<<endl;
			return;
		}
		double r = bef[i]/diag[i-1];
		diag[i] = diag[i] - r*aft[i-1];
		rhsol[i] = rhsol[i] - r*rhsol[i-1];
	}
	rhsol[n-1] = rhsol[n-1]/diag[n-1];
	for(int i=n-2;i>=0;i--)
	{
		rhsol[i] = (rhsol[i] - aft[i]*rhsol[i+1]) / diag[i];
	}
	return;
}

int main()
{
	int kdim = 40;
	double *diag  = new double[kdim];
	double *bef   = new double[kdim];
	double *aft   = new double[kdim];
	double *rhsol = new double[kdim];
	
	int n=6;
	diag[0]  =  1.e0;
	aft[0]   =  0.e0;
	rhsol[0] = 10.e0;
	
	for(int i=1;i<n-1;i++)
	{
		diag[i]  = -2.e0;
		bef[i]	 =  1.e0;
		aft[i]	 =  1.e0;
		rhsol[i] =  0.e0;
	}
	diag[n-1]  =  1.e0;
	bef[n-1]   =  0.e0;
	rhsol[n-1] = 10.e0*n;
	trdiag(kdim,n,bef,diag,aft,rhsol);
	for(int i=0;i<n;i++)
		cout<<i<<"\t"<<rhsol[i]<<endl;
	
	return 0;
}