#include <stdio.h>
#include <stdlib.h>
#include <math.h>
	
void thomas(double a[], double b[], double c[], double r[], int first, int n)
{
/*
Tri-diagonal matrix solver
formated as in fortran source code.  Int(first) has been added so that
one can begin at any starting number (usually 0 or 1) 
Be sure of your precision
implicit double precision (a-h,o-z)
	*/
	int i, im, j, jp, iend;
	iend = first + n - 1;
	c[first] = c[first]/b[first];
	r[first] = r[first]/b[first];
	for (i = (first+1); i <= iend; i++)
	{
		im = i-1;
		c[i] = c[i]/(b[i]-a[i]*c[im]);
		r[i] = (r[i]-a[i]*r[im])/(b[i]-a[i]*c[im]);
	}
	for ( j = iend-1;j >= first; j--)
	{
		jp = j+1;
		r[j] = r[j] - c[j]*r[jp];
	}
}