牛顿迭代法解非线性方程组

#include “stdafx.h“
#include
#include
#define n 101 
#define vt 0.026
#define ni 1e10
#define Na 1e17
#define Nd 2*1e17
#define q 1.6*1e-19
#define E 1e-12
#define h 3*1e-7
using namespace std;
void computeRHS（double x[n]double F[n]）
{int i;
	F[0]=x[0]+vt*log（Na/ni）;
	for（i=1;i<51;i++）
	F[i]=（x[i+1]+x[i-1]-2*x[i]）/pow（h2）-（q/E）*（ni*exp（x[i]/vt）-ni*exp（-x[i]/vt）+Na）;
	for（i=51;i<100;i++）
	F[i]=（x[i+1]+x[i-1]-2*x[i]）/pow（h2）-（q/E）*（ni*exp（x[i]/vt）-ni*exp（-x[i]/vt）-Nd）;
	F[100]=x[100]-vt*log（Nd/ni）;
}
void computeJ（double x[n]double J[n][n]）
{int i;
	for（i=0;i<=100;i++）
		for（int j=0;j<=100;j++）
		J[i][j]=0;
	J[0][0]=1;
	for（i=1;i<=99;i++）
		{
			J[i][i-1]=1/pow（h2）;
			J[i][i]	= -2/pow（h2）-（q/（E*vt））*（ni*exp（x[i]/vt）+ni*exp（-x[i]/vt））;
			J[i][i+1]=1/pow（h2）;
		}
	J[100][100]=1;
}
void inv_J（double J[n][n]double inv[n][n]）
{
	const int n2=n*2;
	double aug[n][n2]L;
	int ijk;
	for （i=0;i	{  for（j=0;j		 aug[i][j]=J[i][j];
	   for（j=n;j		if（j==i+n） aug[i][j]=1;
		else  aug[i][j]=0;
	}
for （i=0;i	{
	  for （k=i+1;k	  {L=-aug[k][i]/aug[i][i];
		for（j=i;j         aug[k][j]=aug[k][j]+L*aug[i][j];
	  }
	}
for （i=n-1;i>0;i--）
	{  
	 for （k=i-1;k>=0;k--）
	  {L=-aug[k][i]/aug[i][i];
		for（j=n2-1;j>=0;j--）
         au