# Using functions instead of lists in physics numerical problems

I'm using the following code for calculation of reflection and transmission coefficients through square well potential:

minE = 5;
maxE = 200;
nE = 200;
energies = Subdivide[minE, maxE, nE - 1];
reflections = ConstantArray[0, nE];
transmissions = ConstantArray[0, nE];

k = Sqrt[2 energies];
Do[
solution =
NDSolve[{psi''[x] == -2 (energies[[i]] - v[x]) psi[x],
psi[xMax] == 1, psi'[xMax] == I*k[[i]]}, psi, {x, -xMax, xMax}];
max = FindMaximum[
Abs[psi[x]] /. solution, {x, -.8 xMax, -.2 xMax}][];
min = FindMinimum[
Abs[psi[x]] /. solution, {x, -.8 xMax, -.2 xMax}][];
reflections[[i]] = ((max - min)/(max + min))^2;,
{i, nE}]

reflFunc = ListInterpolation[reflections, {{minE, maxE}}];
transFunc = ListInterpolation[1 - reflections, {{minE, maxE}}];
Plot[{reflFunc[x], transFunc[x]}, {x, minE, maxE}]


Now, this is kind of problem where I loop over some variable (energy) with manually defined accuracy (by nE and Subdivide). There are two questions:

1. Is there a way to do it with functions instead of loop and lists (how to directly create function I finally get through ListInterpolation) and how to control accuracy/ "energy step"/nE in that case

2. If there is another more elegant way, I'd be most obliged to see it

• As far as I know, for 1D Schrödinger square potential problem, one can find the analytical expression of reflection/transmission as a function of energy. Wolfram can facilitate the derivation work. – Αλέξανδρος Ζεγγ Jul 23 at 0:25

After some reading I'm down to following code, it somehow gets everything into function, but it is surely not most elegant:

ffunc[energy_] := (
k = Sqrt[2 energy];
solution =
NDSolve[{psi''[x] == -2 (energy - v[x]) psi[x], psi[xMax] == 1,
psi'[xMax] == I*k}, psi, {x, -xMax, xMax}];
max = FindMaximum[
Abs[psi[x]] /. solution, {x, -.8 xMax, -.2 xMax}][];
min = FindMinimum[
Abs[psi[x]] /. solution, {x, -.8 xMax, -.2 xMax}][];
((max - min)/(max + min))^2)

v[x_] := If[x > 0 && x < 1, 50, 0];
xMax = 10;
minE = 5;
maxE = 200;
nE = 200;
energies = Subdivide[minE, maxE, nE - 1];
reflections = ffunc /@ energies;
reflFunc = ListInterpolation[reflections, {{minE, maxE}}];
transFunc = ListInterpolation[1 - reflections, {{minE, maxE}}];
ListPlot[reflections]



Anyway, this example evaluates in 9.36s which is a bit longer than original code at 9.02s. Plot with MaxRecursion->0 takes 10.86

I am aware this change does nothing functional, but at least body of Do command is in function now.