How to calculate a gradient at x=0?

Question

I am solving a second order differential equation for n[x]. I solved it with both DSolve and NDSolve. A plot of n as a function of x looks great. Than I want to calculate a gradient of n[x] at a point x=0 and plot it as a function of V. I don't know how to execute this. I calculated gradient at every point and plot it. I get expected shape of the curve and there is also V on x axis, but the values of y axis are to high. Here is the code:

eq = Gx + D[Dif*\!\(
\*SubscriptBox[\(\[PartialD]\), \(x\)]\(n[x]\)\), x] - 
    Kr (n[x] - n0) == 0;
bcs = {n'[d] == 0, n[0] == n0*Exp[(e \[Alpha] V)/(kb T)]};
\[Alpha] = 0.25;
e = 1.602*10^-19;
kb = 1.38*10^-23;
T = 300;
d = 20*10^-6;
Gx = 10^17;
Dif = 4*10^-9;
Kr = 3.1*10^-9;
n0 = 10^21;
sol := {V, 
  D[#, x] &@
   NDSolve[{eq, bcs}, n[x], {x, 0, d}, MaxSteps -> 50000][[1, 1, 2]]}
res = Table[sol, {V, 0, 0.7, 0.05}];
ListPlot[Table[res, {x, 0, d}], Joined -> True, Frame -> True, 
 ImageSize -> 500]

Gx is a complicated function of x because I integrate a list as input. So, my result of DSolve is very big.

I need to stop the calculation of gradient as soon as it becomes zero. When the plot of n'[x] at x=0 cross the x axis (with values of V on it) my calculation should stop and the appropriate result (graph) should be displayed. I don't need the negative values of the derivative, so there is no need for Mathematica for further calculations. What is the best way to execute this?

But I have another problem there because my result is to high and it never changes the sign. In theory it should..

My result should look like:

Because I cannot post the images, you can see the desired graph on the link: https://www.google.si/search?hl=sl&rlz=1C1TEUA_enSI493SI493&q=current+voltage+characteristics+DSSC&bav=on.2,or.r_qf.&biw=1600&bih=785&um=1&ie=UTF-8&tbm=isch&source=og&sa=N&tab=wi&ei=YspCUb3GFYvUsgaTnYGoCA#imgrc=wX5CgAQmm0sdkM%3A%3BJbeOH9F0mOvjNM%3Bhttp%253A%252F%252Fars.els-cdn.com%252Fcontent%252Fimage%252F1-s2.0-S0025540812004047-gr9.jpg%3Bhttp%253A%252F%252Fwww.sciencedirect.com%252Fscience%252Farticle%252Fpii%252FS0025540812004047%3B384%3B314

On x axis is voltage, on y axis is current density, which I should calculate from the gradient at a point x=0: J=(dn[x]/dx)x=0. The current is related to a density gradient, the diffusion coefficient being the proportionality constant. But if I multiply the resulted derivative with Dif, it doesn't change enough. If I multiply the whole equation with e, which gives me correct units, the result is too small.

From the picture above is evident, that my values on y axis are changing to slowly.

Thank you for your answers!

Answer 1

Please try the following

eq = Gx + D[Dif*n'[x], x] -  Kr (n[x] - n0) == 0
bcs = {n'[d] == 0, n[0] == n0*Exp[(e α V)/(kb T)]};
α = 0.25;
e = 1.602*10^-19;
kb = 1.38*10^-23;
T = 300;
d = 20*10^-6;
Gx = 10^17;
Dif = 4*10^-9;
Kr = 3.1*10^-9;
n0 = 10^21;
sol := {V, 
D[#, x] &@NDSolve[{eq, bcs}, n[x], {x, 0, d}, MaxSteps -> 50000][[1, 1, 2]]}
res = Table[sol, {V, 0, 0.7, 0.05}];
Needs["PlotLegends`"]
ListPlot[(TakeWhile[#, #[[2]] > 0&]&/@(Table[#, {x, 0, d, d/100}]&/@Thread[{x, res[[All, 2]]}]))/.{}->Sequence[], Joined -> True, Frame -> True, ImageSize -> 500, PlotLegend -> res[[All, 1]], LegendLabel -> "V", AxesLabel -> {"x", "y"}]