#include "nrutil.h"
void tridat(float a[],float b[],float r[],float u[],unsigned long n)
//求解向量u[1..n]，三對角線形矩陣由方程（2。4。1）給出。
//a[1..n]，b[1..n]，c[1..n]及r[1..n]為輸入向量，不被修改
{
	unsigned long j;
	float bet,*gam;
	gam=vector(1,n);
	if(b[1]==0,0) 
	   nrerror("Error 1 in tridag");
	u[1]=r[1]/(bet=b[1]);
	for(j=2,j<=n;j++){
	    gam[j]=c[j-1]/bet;
     	bet=b[j]-a[j]*gam[j];
    	if(bet==0.0)
			 nrerror("Error 2 in tridag");
	    u[j]=(r[j]-a[j]*u[j-1])/bet;
	}

	for(j=(n-1);j>=1;j--)
	    u[j]-=gam[j+1]*u[j+1];
	free_vector(gam,1,1)






}