# How to calculate a gradient at x=0?

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 == 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:

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!

• Why do you say it is too high? – wxffles Mar 13 '13 at 21:45
• Your code seems right. If your result for the derivative are much higher than you expected, perhaps you could check those constants. – Dr. belisarius Mar 13 '13 at 23:13
• You can see analytically from the DSolve solution that this is correct. Why do you think that this is the wrong answer? – Jonathan Shock Mar 14 '13 at 6:29
• I want to observe a current density. It is the gradient that I am calculating. But I should multiply my result by a proportionality factor Dif, which changes the result but not enough. The current density cannot have so high values. It is true, that also the units are not correct. I should multiply the result for gradient by elementary charge to get the current. But even I do this, the values are not ok. I must also find a way to stop the calculation, when my gradient (current density) is zero. I think If loop will do it. But first, I want to be sure, my calculation is correct.. – Luka Mar 14 '13 at 6:57
• You can use the command AspectRatio if you just want to change the scale between the axes, for more see here – Spawn1701D Mar 14 '13 at 19:09

## 1 Answer

Please try the following

eq = Gx + D[Dif*n'[x], x] -  Kr (n[x] - n0) == 0
bcs = {n'[d] == 0, n == 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[#, #[] > 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"}]
`