пробирает страх

[ciaonataha]

Есть функция equilibre_Fe_Cr_C, которая зависит от вектора Y. Этот вектор я определяю в самой функции для pot_chemical1[i] ,pot_chemical2[i] + определяю в main.cpp....
Значения в векторе Y будут изменяться после каждой итерации за счет метода Ньютона Рафсона. Так вот. Меня пробирает страх от того, что не будет пересчитиваться pot_chemical1, pot_chemical2 для новых значении вектора. Сразу хочу сказать, что функцию, где я записиваю нельзя сделать зависимой от вектора.Жду с нетерпением предложений

void interface::equilibre_Fe_Cr_C(int n, std::vector <double> &Y, std::vector <double> &F)
{
Y[0] = ph1->Thermo->y1[0]; //printf("%lf\n",Y[0]);
Y[1] = ph1->Thermo->y1[1]; //printf("%lf\n",Y[1]);
Y[2] = ph1->Thermo->y2[0]; //printf("%lf\n",Y[2]);
Y[3] = ph1->Thermo->y2[1]; //printf("%lf\n",Y[3]);
Y[4] = ph2->Thermo->y1[0]; //printf("%lf\n",Y[4]);
Y[5] = ph2->Thermo->y1[1]; //printf("%lf\n",Y[5]);
Y[6] = ph2->Thermo->y2[0]; //printf("%lf\n",Y[6]);
Y[7] = ph2->Thermo->y2[1]; //printf("%lf\n",Y[7]);
Y[8] = X_nf;  //printf("%lf\n",Y[8]);
//printf("%lf\n",X_C);
//printf("%lf\n",X_Cr);

 F[0] = ph1->Thermo->pot_chemical1[0] - ph2->Thermo->pot_chemical1[0];// printf("%lf\n",F[0]);
 F[1] = ph1->Thermo->pot_chemical1[1] - ph2->Thermo->pot_chemical1[1];// printf("%lf\n",F[1]);
 F[2] = ph1->Thermo->pot_chemical2[0] - ph2->Thermo->pot_chemical2[0];//printf("%lf\n",F[2]);
 F[3] = Y[0] + Y[1] - 1; //printf("%i\n",F[3]);
 F[4] = Y[4] + Y[5] - 1; //printf("%i\n",F[4]);
 F[5] = Y[2] + Y[3] - 1; //printf("%i\n",F[5]);
 F[6] = Y[6] + Y[7] - 1; //printf("%i\n",F[6]);
 F[7] = (ph1->Thermo->m[0] * Y[1])/(ph1->Thermo->m[0] + ph1->Thermo->m[1] - ph1->Thermo->m[1] * Y[3]) * Y[8] + (ph2->Thermo->m[0] * Y[5])/(ph2->Thermo->m[0] + ph2->Thermo->m[1] - ph2->Thermo->m[1] * Y[7]) * (1 - Y[8]) - X_Cr; //printf("%lf\n",F[7]);
 F[8] = (ph1->Thermo->m[1] * Y[2])/(ph1->Thermo->m[0] + ph1->Thermo->m[1] - ph1->Thermo->m[1] * Y[3]) * Y[8]  + (ph2->Thermo->m[1] * Y[6])/(ph2->Thermo->m[0] + ph2->Thermo->m[1] - ph2->Thermo->m[1] * Y[7]) * (1 - Y[8]) - X_C; //printf("%lf\n",F[8]);
}

void Gibbs::chemical_potential()
{
 int i, j;
 int i1 = dGm1.size();
 int j1 = dGm2.size();
 double dummy1 = 0., dummy2 = 0.;

 for (i = 0; i < i1; i++){
  dummy1 += y1[i] * dGm1[i];}
 for (i = 0; i < i1; i++){
  pot_chemical1[i] = Gm + 1/m[0]*(dGm1[i] - dummy1);} 

 for (j = 0; j < j1; j++){
  dummy2 += y2[j] * dGm2[j];}
 for (j = 0; j < j1; j++){
  pot_chemical2[j] = Gm + 1/m[1]*(dGm2[j] - dummy2);}
}

[("")(Э_Є)("")]

А меня пробирает страх от "Ньютона Рафсона" и мурашки по телу от листингов. Проблема то в чем? То есть оно не работает? А можно увидеть формулу ?

[ciaonataha]

проблемма в том , что void chemical_potential() не зависит на прямую от вектора Y. Она зависит от состовляющих вектора Y. Вот я и не знаю, будет ли пересчитиваться функция void chemical_potential() после каждой итерации. Пока что я не дошла до этой ошибки. Просто, предчувствие подсказывает что здесь что-то не так. Спрашываю совета у публики.