#pragma rtGlobals=3		// Use modern global access method and strict wave access.

Function EquSystem(pw, xx, yw, dydx)
		Wave pw // constants:pw[1], ..., pw[5]
		Variable xx
		Wave yw // yw[0]=n, yw[1]=p, yw[2]=jn, yw[3]= jp, yw[4]=p (potential), yw[5]=f (electric field)
		Wave dydx

		dydx[0]=pw[0]*yw[2]-pw[2]*yw[0]*yw[5]
		dydx[1]=pw[2]*yw[1]*yw[5]-pw[1]*yw[3]
		
		dydx[2]=-(pw[3]-pw[4]*yw[0]*yw[1])
		dydx[3]=pw[3]-pw[4]*yw[0]*yw[1]
		
		dydx[4]=yw[5]
		dydx[5]=pw[5]*(yw[1]-yw[0])
		
		
		
		//pw[0]=(1/(q*Dn))=(1/(mn*kb*T))  -  mn=2.5*(10^(-7)) (m^2)/(V*s), kb=1.38*10^(-23), T=300 K  -   9.66184*10^26
		//pw[1]=(1/(q*Dp))=(1/(mp*kb*T))  -  mp=3*(10^(-8))  -  8.05153*10^27
		//pw[2]=(mn/Dn)=(1/kb*T)  -   2.41546*10^20
		//pw[3]=q*G  -  q=1.6*(10^(-19)), , G=2.7*10^27  -  4.32*10^8
		//pw[4]=((q^2)*(mn+mp)/e)   -   7.168*10^(-45)
		//pw[5]=(-(q/e0))  -  e0=3*10^(-11)  -  5.33*10(-31)


		return 0
End


Function PlotJV()
Variable  i=0, dim=11  // dim - pionts quantity for viltage

//initial conditions for:
Variable n1=2.5e25, n2=4.11e8			//  electron density
//Variable p1, p2		//  hole density
Variable jn1=0.18e-2, jn2=2.4e-2	//  current density of electrons
Variable jp1=0.62e-2, jp2=0	//  current density of holes
Variable pot1, pot2 	//  potential 
Variable f1=2.7e7, f2=0            // electric field

// n1=Nc 
// n2=Nc*exp(-Egap*q/kb*T)  -  Egap=1eV=q
// p1=Nv*exp(-Egap*q/kb*T)
// p2=Nv

make/O/D/N=(dim) v	// voltage

// vbi - built-in-potential  voltage = band gap energy
Variable vbi=1
v[0]=0.1


	for(i=0;i<dim;i+=1)
		v[i]=(i+1)*0.1
		pot1=((vbi-v[i])/2)
		pot2=((v[i]-vbi)/2)

	endfor	
	
		Make/D/O/N=(90,6) Equat
		SetDimLabel 1,0,n,Equat
		SetDimLabel 1,1,p,Equat
		SetDimLabel 1,2,jn,Equat
		SetDimLabel 1,3,jp,Equat
		SetDimLabel 1,4,pot,Equat
		SetDimLabel 1,5,f,Equat

		Equat[0][%n]=n1
		Equat[0][%p]=n1
		Equat[0][%jn]=jn1
		Equat[0][%jp]=jp1
		Equat[0][%pot]=pot1
		Equat[0][%f]=f1
			
	
		//Make/D/O wave0={1e24, 1e28, 1e6, 1e6, 1e-33, 1e-34}
		//Make/D/O wave0={1e-3, 1e-3, 1e-3, 1e-3, 1e-3, 1e-3}
		Make/D/O wave1={n2, n2, jn2, jp2, pot2, f2} // conditions
		Make/D/O PP={9.66184e26,  8.05153e27,  2.41546e20,  4.32e8,  7.168e-45, 5.33e-31}
		IntegrateODE/M=2/STOP={wave1,0} EquSystem, PP, Equat		
		
		//S=wave0/F=(1+4)/E=1e-3/
		//CVOP={0,0,1,1e-7}/
		
		Duplicate/O Equat $"Equat_V_"+num2str(v[i])
		
		//Print i
		//Display v

//	Make/D/O/N=10 current
//	SetFormula current, "x"
//	Make/D/O/N=10 voltage
//	SetFormula voltage, "y"
//	Display current vs voltage	

	Display Equat[][%jn]
	Print V_ODEStepSize

return 0

end