#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