pls help me i have a problem as listed below:
i had used lud, lubksb and mprove routines in NR C++, but i can't get a physical result:
every element in the matrices is in equation forms or zero value. Like A is
dxdy[1] = al+3.0*beta*gx[1]*gx[1]+beta*gy[1]*gy[1];
dxdy[2] = j1/2.0+j2/8.0;
dxdy[3] = j1/2.0+j2/8.0;
dxdy[4] = 2.0*beta*gy[1]*gx[1];
dxdy[5] = 0;
dxdy[6] = 0;

dxdy[7] = j1/2.0+j2/8.0;
dxdy[8] = al+3.0*beta*gx[2]*gx[2]+beta*gy[2]*gy[2];
dxdy[9] = j1/2.0+j2/8.0;
dxdy[10] = 0;
dxdy[11] = 2.0*beta*gy[2]*gx[2];
dxdy[12] = 0;

dxdy[13] = j1/2.0+j2/8.0;
dxdy[14] = j1/2.0+j2/8.0;
dxdy[15] = al+3.0*beta*gx[3]*gx[3]+beta*gy[3]*gy[3];
dxdy[16] = 0;
dxdy[17] = 0;
dxdy[18] = 2.0*beta*gy[3]*gx[3];

dxdy[19] = 2.0*beta*gy[1]*gx[1];
dxdy[20] = 0;
dxdy[21] = 0;
dxdy[22] = al+3.0*beta*gy[1]*gy[1]+beta*gx[1]*gx[1];
dxdy[23] = j1/2.0+j2/8.0;
dxdy[24] = j1/2.0+j2/8.0;

dxdy[25] = 0;
dxdy[26] = 2.0*beta*gy[2]*gx[2];
dxdy[27] = 0;
dxdy[28] = j1/2.0+j2/8.0;
dxdy[29] = al+3.0*beta*gy[2]*gy[2]+beta*gx[2]*gx[2];
dxdy[30] = j1/2.0+j2/8.0;

dxdy[31] = 0;
dxdy[32] = 0;
dxdy[33] = 2.0*beta*gy[3]*gx[3];
dxdy[34] = j1/2.0+j2/8.0;
dxdy[35] = j1/2.0+j2/8.0;
dxdy[36] = al+3.0*beta*gy[3]*gy[3]+beta*gx[3]*gx[3];



//declaration of 6X6 matrix.

const DP a_d[N*N]=

{

dxdy[1], dxdy[2], dxdy[3], dxdy[4], dxdy[5], dxdy[6],

dxdy[7], dxdy[8], dxdy[9], dxdy[10], dxdy[11], dxdy[12],

dxdy[13], dxdy[14], dxdy[15], dxdy[16], dxdy[17], dxdy[18],

dxdy[19], dxdy[20], dxdy[21], dxdy[22], dxdy[23], dxdy[24],

dxdy[25], dxdy[26], dxdy[27], dxdy[28], dxdy[29], dxdy[30],

dxdy[31], dxdy[32], dxdy[33], dxdy[34], dxdy[35], dxdy[36]

};

//And the B is

de[1] = - ((al*gx[1])+(beta*gx[1]*(gx[1]*gx[1]+gy[1]*gy[1]))+(j1/4.0*(gx[2]+gx[3]))+(j2/8.0*(gx[2]+gx[3])));
de[2] = - ((al*gx[2])+(beta*gx[2]*(gx[2]*gx[2]+gy[2]*gy[2]))+(j1/4.0*(gx[1]+gx[3]))+(j2/8.0*(gx[1]+gx[3])));
de[3] = - ((al*gx[3])+(beta*gx[3]*(gx[3]*gx[3]+gy[3]*gy[3]))+(j1/4.0*(gx[1]+gx[2]))+(j2/8.0*(gx[1]+gx[2])));
de[4] = (-1.0*E) - ((al*gy[1])+(beta*gy[1]*(gx[1]*gx[1]+gy[1]*gy[1]))+(j1/4.0*(gy[2]+gy[3]))+(j2/8.0*(gy[2]+gy[3])));
de[5] = (-1.0*E) - ((al*gy[2])+(beta*gy[2]*(gx[2]*gx[2]+gy[2]*gy[2]))+(j1/4.0*(gy[1]+gy[3]))+(j2/8.0*(gy[1]+gy[3])));
de[6] = (-1.0*E) - ((al*gy[3])+(beta*gy[3]*(gx[3]*gx[3]+gy[3]*gy[3]))+(j1/4.0*(gy[1]+gy[2]))+(j2/8.0*(gy[1]+gy[2])));


const DP b_d[N]={ de[1], de[2], de[3], de[4], de[5], de[6]};



cout<<endl;


DP d;
Vec_INT indx(N);
Vec_DP x(N),b(b_d,N);
Mat_DP aa(N,N),a(a_d,N,N);

for (i=0;i<N;i++) {
x[i]=b[i];
for (j=0;j<N;j++)
aa[i][j]=a[i][j];
}
NR::ludcmp(aa,indx,d);
NR::lubksb(aa,indx,x);
cout << endl << "Solution vector for the equations:" << endl;
for (i=0;i<N;i++) cout << setw(12) << x[i];
cout << endl;



for(i=1; i<=10; i++)
{
NR::mprove(a,aa,indx,b,x);
}

cout << endl << "Solution vector recovered by mprove:" << endl;
for (i=0;i<N;i++) cout << setw(12) << x[i];
cout << endl;


there is one physical argument, if i let the value of E in the B matrix to be zero, then all the results from x[i] should be zero also, but this is not the case when i use the NR subroutines. May i know the reason? i had tried to use the mathematica to solve the same problem and it gives the sensible result.
i still need to use c++ because it is a more suitable tool to solve the problems i involve now. Waiting for favourable reply. :)

The matrices are from a function, let say f(x)
{Hessian of f(x)}.{X}={B - first derivative of f(x)}
so when B=0, then the X will be resulted as zero. but when i let the B=0, there are non-zero values in X[i]. i use a unix system, something wrong?
hopefully someone can help me.

:)

after i check carefully, then it works ok now, but still have some problem.