/** Fills solution into x. Warning: will modify c and d! **/
void TridiagonalSolve(const double *a, const double *b, double *c,
			    double *d, double *x, unsigned int n)
{
	int i;

	// Modify the coefficients.
	c[0] = c[0]/b[0];			       //Division by zero risk.
	d[0] = d[0]/b[0];
	double id;
	for(i = 1; i != n; i++){
		id = 1.0/(b[i] - c[i - 1]*a[i]);	//Division by zero risk.
		c[i] = c[i]*id;			 //Last value calculated is redundant.
		d[i] = (d[i] - a[i]*d[i - 1])*id;
	}

	//Now back substitute.
	x[n - 1] = d[n - 1];
	for(i = n - 2; i != -1; i--)
		x[i] = d[i] - c[i]*x[i + 1];
}