# Defined function keeps the variable

I'm a beginner at Mathematica, and I have a problem finding the roots of a function for a specific variable (using Jens' findAllRoots function).

Due to the nature of my problem, I need to compute these zeros inside a table on the fly. I decided to define a function to calculate them, which depends on a variable. Nevertheless, after some iterations, the function keeps information from the previous value given to the variable. What is going wrong? Thanks for your help.

Here is a sample code:

ClearAll[V, g, tem, lim1, fun]

V = 100;
g = 3;
lim1 = 20;

fun[x_,tem_] = D[
TT[ x , g ] * fte[ x , +V/2, tem]*( 1 - fte[ x , +V/2, tem]),
x ]; (* TT is a transmission function and fte is the Fermi distribution  *)

zeros[tem_] := findAllRoots[
fun[x , tem ], { x , -V/2 - 20*tem, V/2 + 20*tem},
"ShowPlot" -> False]


The derivative (function fun) should be computed only once, keeping x and tem as parameters. The function zeros must be computed every time it is called. My first problem is that zeros[tem] does not work, it does not provide any number. To solve it, I have tried using:

zeros1 = Function[tem,
findAllRoots[fun[ x , tem ], { x , -V/2 - 20*tem, V/2 + 20*tem},
"ShowPlot" -> False]
]

zeros2 =   (findAllRoots[fun[ x , #], {x, -V/2 - 20*#, V/2 + 20*#}, "ShowPlot" -> False]) &


But after some iterations the functions maintain the previous given value

zeros1[45] = {-10.8508, 4.43734*10^-31, 50.562}

zeros1[50] = {-10.8508, 4.43734*10^-31, 50.562}[50]

What is the best form to define this zeros function?

(* The definition of the other functions *)

TT[x_, g_] := e^2/(e^2 + g^2)

fte[x_, m_, t_] := 1/(Exp[(e - m)/t] + 1)

SyntaxInformation[findAllRoots] = {"LocalVariables" -> {"Plot", {2, 2}},
"ArgumentsPattern" -> {_, _, OptionsPattern[]}};
SetAttributes[findAllRoots, HoldAll];

Options[findAllRoots] = Join[{"ShowPlot" -> False, PlotRange -> All},
FilterRules[Options[Plot], Except[PlotRange]]];

findAllRoots[fn_, {l_, lmin_, lmax_}, opts : OptionsPattern[]] :=
Module[{pl, p, x, localFunction, brackets}, localFunction = ReleaseHold[Hold[fn] /. HoldPattern[l] :> x];

If[lmin != lmax, pl = Plot[localFunction, {x, lmin, lmax},
Evaluate@FilterRules[Join[{opts}, Options[findAllRoots]], Options[Plot]]];

p = Cases[pl, Line[{x__}] :> x, Infinity];

If[OptionValue["ShowPlot"], Print[Show[pl, PlotLabel -> "Finding roots for this function", ImageSize -> 200, BaseStyle -> {FontSize -> 8}]]], p = {}];

brackets = Map[First, Select[(*This Split trick pretends that two points on the curve are "equal" if the function values have _opposite _ sign.Pairs of such sign-changes form the brackets for the subsequent FindRoot*)

Split[p, Sign[Last[#2]] == -Sign[Last[#1]] &], Length[#1] == 2 &], {2}];

x /. Apply[FindRoot[localFunction == 0, {x, ##1}] &, brackets, {1}] /. x -> {}]

• Please provide the definitions for all functions (e.g., TT, fte, findAllRoots) Commented Sep 22, 2020 at 21:06
• I assume that the definitions for TT and fte should have x instead of e on the right hand side? Also, why can't you just use Solve? Commented Sep 22, 2020 at 21:45
• Yes, it was a mistake while I was writing the post. If I am right, I tried using solve as first option, but the answer didn't provide all the possible solutions. Commented Sep 23, 2020 at 15:36

Clear["Global*"]

TT[x_, g_] := x^2/(x^2 + g^2)

fte[x_, m_, t_] := 1/(Exp[(x - m)/t] + 1)

V = 100;
g = 3;
lim1 = 20;

fun[x_, tem_] =
D[TT[x, g]*fte[x, +V/2, tem]*(1 - fte[x, +V/2, tem]), x] // FullSimplify
(*TT is a transmission function and fte is the Fermi distribution*)

(* -((x Sech[(-50 + x)/(
2 tem)]^2 (-18 tem + x (9 + x^2) Tanh[(-50 + x)/(2 tem)]))/(
4 tem (9 + x^2)^2)) *)


Using NSolve to find the roots

zeros1[tem_?NumericQ] :=
NSolve[{fun[x, tem] == 0, -V/2 - 20*tem <= x <= V/2 + 20*tem}, x]

zeros1[45]

(* {{x -> -10.8508}, {x -> 0.}, {x -> 50.562}} *)


x == 0 is always a root

fun[0, tem]

(* 0 *)

zeros1[50]

(* {{x -> -11.5445}, {x -> 0.}, {x -> 50.6887}} *)


To plot zeros1

data = Transpose[
Table[Thread[{tem, x /. zeros1[tem]}], {tem , 5/2, 50, 5/2}]];

MinMax[#[[All, 2]]] & /@ data

(* {{-11.5445, -2.73266}, {0., 0.}, {50.0018, 50.6887}} *)

ListLinePlot[data,
Frame -> True,
FrameLabel -> (Style[#, 14, Bold] & /@ {tem, x}),
PlotLegends -> Automatic,
PlotLabel -> Style["fun[x, tem] \[Equal] 0", 14, Bold]]


Comparing with ContourPlot

ContourPlot[fun[x, tem] == 0, {tem, 5/2, 50}, {x, -15, 54},
PlotRange -> {{-1, 51}, {-15, 54}},
WorkingPrecision -> 50,
PlotPoints -> 50,
MaxRecursion -> 5,
FrameLabel -> (Style[#, 14, Bold] & /@ {tem, x}),
PlotLabel -> Style["fun[x, tem] \[Equal] 0", 14, Bold],
AspectRatio -> 1/GoldenRatio]
`

• Really appreciate it! This works nicely. Commented Sep 25, 2020 at 18:43