I have some code which numerically solves a differential equation for a given value of the parameter $M0$ and for some initial conditions. The initial conditions will remain fixed throughout the procedure and we want to vary the value of the parameter $M0$.
The code is given below
eqnofmotion := 1/r^3 D[r^3 D[f[r], r], r] + M0^2/(r^2 + L^2)^2 f[r]
L := 1
M0 := 2 Sqrt[2]
sltn = NDSolve[{eqnofmotion == 0, f[10^-5] == 1, f'[10^-5] == 0},
f, {r, 10^-5, 10^2}];
Plot[ f[r] /. sltn, {r, 10^-5, 10^2}, PlotRange -> {-1, 1},
PlotStyle -> {Thick, Black},
BaseStyle -> {18, FontFamily -> "Times New Roman"},
AxesLabel -> {"\[Rho]", "\[Delta]L(\[Rho])"},
PlotLegends -> Automatic]
f[r] /. sltn /. r :> 10^2
After this, the following procedure takes place. Define as a function the solution above, and its derivative and evaluate them at $100$.
function0[r_] := f[r] /. sltn
derivative0[r_] := function0'[r]
function0[r] /. r :> 100
derivative0[r] /. r :> 100
From the above values, we can calculate the operator and the source of the state. The code is given below
scalarasymptote[r_] := operator/r^3 + source r
(*Define and solve the system to determine operator and source*)
\
Solve[{scalarasymptote[100] == 0.0000999607,
scalarasymptote'[100] == - 1.99889 10^-6}, {operator, source}]
And now that we have the values of the operator and the source, we can calculate the square of the coupling, which is the $a$ in the following code,
(*Now use that \[ScriptCapitalJ] = g^2/\[CapitalLambda]^2 \
\[ScriptCapitalO]*)
(*\[CapitalLambda] is the UV cutoff.*)
\[CapitalLambda] := 10^2
NSolve[ 2.49983 10^-7 == 74.9624 a/\[CapitalLambda]^2, a]
The goal is to repeat all of the above steps for many values of $M0$, that are smaller than $2 \sqrt{2}$. Create a list of the square of all the values of $M0$ and the corresponding $a$ values and plot the data.
The way I do it, is to copy and paste the above piece of code, and run it each time for a different value of $M0$, for which I use another symbol to save it, find $a$ each time, until I reach the value $M0=0$ and then plot the data.
The thing I would like to know is if and how I can automate the above procedure or at least some steps of it, in order to get more points a bit more easily and within a reasonable amount of time. Copying and pasting and running the whole thing each time takes quite some time and the goal is to do this for excited states as well, start with higher $M0$ and go all the way down to $0$.
Thanks in advance.