Hi,

Following is the segment of my main C code:
poly_fit () is my own function.
I am using Numerical Recipe subroutine gaussj.c and function dmatrix() from nrutil.c

Basically I want to solve linear equations.
I am fitting a 4th degree polynomial fucntion to some data points.

I get matrix equation: A*B=C
where A= coef_matrix in my code
B= set of 5 coefficients to be found e.g. {a0,a1,a2,a3,a4} = sol_col in my code
C= data points obtained from array old_est[] as argument to my function = sol_col in my code
(sol_col is replaced by the solution from gaussj.c )

N = 5 (N is global variable, there are other glbal bariables also in this segment).

My problems:
1. Is use of dmatrix() okay?
2. I get the output which contains very large unexpected numbers alternately i.e. 1 st no ok, 2nd no bad, 3 rd no ok,....so on.

I AM SURE THAT PROBLEM IS WITHIN MY poly_fit() FUNCTION REGARDING dmatrix() and gaussj.c

Any help will be appreciated.


Code:
void poly_fit(double ip[], double old_est[], double err[])
 {
	double **coef_matrix, **sol_col;        
       
	coef_matrix=dmatrix(1, N, 1, N);
	sol_col=dmatrix(1, N, 1, 1);

       
    int i,j,n,r,c;
    
    nst_neg=-1*nst;

//-------------resetting to zero ------------------------
    for(r=1;r<=N;r++) 
	{
	   for(c=1;c<=N;c++)
		{
                   coef_matrix[r][c]=0.0;	
		}
	}
//-------------------------------------------------------------

    for(j=nst_neg;j<nst+1;j++)  // j varies from -10 to +10
  
     {
         for(r=1;r<=N;r++)
	   {
	      for(c=1;c<=N;c++)
		{
		   coef_matrix[r][c] = coef_matrix[r][c] + pow( j,(r+c-2)*1.0 );
		}
	   }    
     } 


    for(n=nst+1;n<nend+2;n++)   // n varies from 10 to 290
      {

	    for(r=1;r<=N;r++)
		{
		     c=1;
		    sol_col[r][c]=0;
		}   
          
            for(j=nst_neg;j<nst+1;j++)       // j varies from -10 to +10
             {
              	for(r=1;r<=N;r++)
	   	  {
		    c=1;

 		    sol_col[r][c]= sol_col[r][c] + (old_est[n+j]*pow(j,(r-1)*1.0));
			
		  }    
             }
             
	
//---------------------------------------------------------

      		gaussj(coef_matrix, N, sol_col, 1);      
 /* numerical recipe for solving simultaneous equations */
		
//-----------------------------------------------------------

	     err[n] =sol_col[1][1]- ip[n] ;			
     
       }

             free_dmatrix(coef_matrix, 1, N, 1, N);           
             free_dmatrix(sol_col, 1, N, 1, 1);		
           
    return ;

 }
thanks