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



Function Test(pw, xx, yw, dydx)
		Wave pw // pw[0]=k1, pw[1]=k2, pw[2]=k3, pw[3]=k4
		Variable xx
		Wave yw // yw[0]=n, yw[1]=j, yw[2]=pot, yw[3]=f
		Wave dydx

		dydx[0]=pw[0]*yw[1]-pw[1]*yw[0]*yw[3]
		dydx[1]=pw[2]
		dydx[2]=yw[3]
		dydx[3]=pw[3]*yw[0]
		
		return 0
End


Function GraphJV()
Variable i=0, dim=11, vbi=1
Variable n1=2.5, n2=0.5, j1=0.13, j2=0.01, f1=5, f2=0, pot1, pot2
make/O/D/N=(dim) v


	v[0]=0.1
	for(i=0;i<dim;i+=1)
		//v[i]=0.1*(i+1)
		v[i]=(i+1)*0.1
		pot1=(vbi-v[i])/2
		pot2=(v[i]-vbi)/2
		
	endfor
		
		Make/D/O/N=(90,4) Condit
		SetDimLabel 1,0,n,Condit
		SetDimLabel 1,1,j,Condit
		SetDimLabel 1,2,pot,Condit
		SetDimLabel 1,3,f,Condit
		
		Condit[0][%n]=n1 
		Condit[0][%j]=n2
		Condit[0][%pot]=pot1
		Condit[0][%f]=f1
		//Make/D/O wave0={.1,.01,.1*10^(-5),.1*10^(-3)}
		Make/D/O PP={96,0.024,7.168e-5, 5.33e-3}
		IntegrateODE Test, PP, Condit

		Display Condit[][%j]

return 0

end