μ=0.7;
d=0.01;
sol=NDSolve[{μ ρ'[x]Ex[x]-d ρ''[x]==0,Ex'[x]==0.001ρ[x],ρ[0]==3,ρ[1]==1,Ex[1]==-0.051},{ρ[x],Ex[x]},{x,0,1}];
Plot[ρ[x]/.sol,{x,0,1},PlotRange->All]
Plot[Ex[x]/.sol,{x,0,1},PlotRange->All]

μe=0.7;
μh=0.7;
d=0.01/0.7;
g[x_]=0.05 Exp[-(1/2)(((((x-0.5)/0.1)))^2)];
r=0.02;
sol=NDSolve[{μh p'[x]Ex[x]-μh d p''[x]==-r n[x]p[x]+g[x],-μe n'[x]Ex[x]-μe d n''[x]==-r n[x]p[x]+g[x],Ex'[x]==0.001(p[x]-n[x]),p[1]==0.3,p[0]==0.1,n[0]==0.3,n[1]==0.1,Ex[1]==0.051},{n[x],p[x],Ex[x]},{x,0,1}];
Plot[g[x],{x,0,1},PlotRange->All]
Plot[{(p[x]/.sol),(n[x]/.sol)},{x,0,1},PlotRange->All]
Plot[Ex[x]/.sol,{x,0,1},PlotRange->All]
charge[x_]:=(p[x]-n[x])/.sol[[1]]
Plot[charge[x],{x,0,1}]
NIntegrate[charge[x],{x,0,1}]

μe=1.4;
μh=0.7;
d=0.01/0.7;
g[x_]=0.0 Exp[-(1/2)(((((x-0.5)/0.1)))^2)];
r=0.00002;
data1={};
Do[vol=0.1i/30;
sol=NDSolve[{-μh p'[x]v'[x]-μh p[x]v''[x]-μh d p''[x]==-r n[x]p[x]+g[x],μe n'[x]v'[x]+μe n[x]v''[x]-μe d n''[x]==-r n[x]p[x]+g[x],v''[x]==-0.001(p[x]-n[x]),p[1]==0.3,p[0]==0.1,n[0]==0.3,n[1]==0.1,v[0]==vol,v[1]==0.0},{n[x],p[x],v[x]},{x,0,1}];
n0[x_]:=n[x]/.sol[[1]];
p0[x_]:=p[x]/.sol[[1]];
v0[x_]:=v[x]/.sol[[1]];
(* Plot[g[x],{x,0,1},PlotRange->All]
Plot[{n0[x],p0[x]},{x,0,1},PlotRange->All]
Plot[v0[x],{x,0,1},PlotRange->All] *)
charge[x_]:=(p0[x]-n0[x]);
(* Plot[charge[x],{x,0,1}]
NIntegrate[charge[x],{x,0,1}] *)
currh0=(-μh p0[x]D[v0[x],x]-μh d D[p0[x],x]-(μe n0[x]D[v0[x],x]-μe d D[n0[x],x]))/.{x->0};
curr1=(-μh p0[x]D[v0[x],x]-μh d D[p0[x],x]-(μe n0[x]D[v0[x],x]-μe d D[n0[x],x]))/.{x->1};
data1=Append[data1,{vol,curr1}];
Print[i];
,{i,0,30}]

ListPlot[data1]

μe=1.4;
μh=0.7;
d=0.01/0.7;
g[x_]=0.1 Exp[-(1/2)(((((x-0.5)/0.03)))^2)];
r=0.02;
vol=0.1;
sol=NDSolve[{-μh p'[x]v'[x]-μh p[x]v''[x]-μh d p''[x]==-r n[x]p[x]+g[x],μe n'[x]v'[x]+μe n[x]v''[x]-μe d n''[x]==-r n[x]p[x]+g[x],v''[x]==-0.001(p[x]-n[x]),p[1]==0.1,p[0]==0.1,n[0]==0.1,n[1]==0.1,v[0]==vol,v[1]==0.0},{n[x],p[x],v[x]},{x,0,1}];
n0[x_]:=n[x]/.sol[[1]];
p0[x_]:=p[x]/.sol[[1]];
v0[x_]:=v[x]/.sol[[1]];
charge[x_]:=(p0[x]-n0[x]);
Plot[g[x],{x,0,1},PlotRange->All]
Plot[{n0[x],p0[x]},{x,0,1},PlotRange->All]
Plot[v0[x],{x,0.4,0.6},PlotRange->All]
Plot[charge[x],{x,0,1}]
El=-D[v0[x],x];
Plot[El,{x,0,1}]
NIntegrate[charge[x],{x,0,1}]
curr0=(-μh p0[x]D[v0[x],x]-μh d D[p0[x],x]-(μe n0[x]D[v0[x],x]-μe d D[n0[x],x]))/.{x->0};
curr1=(-μh p0[x]D[v0[x],x]-μh d D[p0[x],x]-(μe n0[x]D[v0[x],x]-μe d D[n0[x],x]))/.{x->1};

El[x]/.{x->0.5}