Defined function keeps the variable

Question

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 -> {}]

Answer 1

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]