void tridiag(double *a, double *b, double *c, double *r, double *u, int n)
{
	int j;
	double bet, gam[n-1];
	/* a is AL, b is AM, c is AR */
	/* r is rhs, u is output, n is dimension */
	
	if(b[0]==0) printf("tridag: rewrite equations");
	bet=b[0];
	u[0]=r[0]/bet;
	for(j=1;j<n;j++){
		gam[j-1]=c[j-1]/bet;
		bet=b[j]-a[j]*gam[j-1];
		if(bet==0) printf("tridag failed");
		u[j]=(r[j]-a[j]*u[j-1])/bet;
	}
	for(j=n-2;j>-1;j--){
		u[j]=u[j]-gam[j]*u[j+1];
	}
	return;
}