Writing NDSolve-like functions

Question

I've been using Mathematica for a while now, and there are multiple common tasks that I end up doing. I'd like to learn how to convert them into functions so that I can store them in some file and load them at startup.

Most of these common tasks involve taking a set of equations and doing something with them. Thing is, I don't want to fix the variables the equations are written in, and would like behavior similar to NDSolve where one can specify the functions as well as which variable is which.

For example, take this task:

Pd = ParametricNDSolve[{v'[t] == (v[t] - v[t]^3/3 - ω[t])/0.2, 
       ω'[t] == (1 + v[t])/0.2, 
       ω[0] == b, v[0] == a}, {v[t], ω[t]}, 
        {t, 0, 20}, {a, b}];
PVec[x_, y_] := {D[(v[t] /. Pd) [x, y], t] /. t -> 0, 
                 D[(ω[t] /. Pd) [x, y], t] /. t -> 0};
VectorPlot[PVec[x, y], {x, -2, 1}, {y, -1.2, 0.5}]

Here, I'm parametrically solving the pair of equations

$$ \begin{align} v'(t)&=\frac{-\frac{1}{3} v(t)^3+v(t)-\omega (t)}{0.2}\\ \omega'(t)&=\frac{v(t)+1}{0.2} \end{align} $$

with the parameters being the initial state of $v$ and $\omega$. Now, I take the vector field of their derivatives at $t=0$, and vector plot it. Basically, I want to get a quick map of how the initial trajectories ($(v'(t),\omega'(t))$) change for different initial conditions without having to open EquationTrekker and do things manually.

---

So now, I want to take this task and write it as a function akin to TrajectoryMap[{eqn1,eqn2},{{var1,varmin,varmax},{var2,varmin,varmax}},{{timevar,timepoint},tmin,tmax,tpoint}] where this particular code would run as

TrajectoryMap[{v'[t] == (v[t] - v[t]^3/3 - ω[t])/0.2, 
    ω'[t] == (1 + v[t])/0.2}, {{v[t],-2,1}, {ω[t],-1.2,0.5}}, {{t,0},0,20}]

(or something similar)

Here, timepoint is the value of t for which the derivative of the interpolating function is calculated.

I can tell that this will probably need some combination of Holds and Evaluates, but I can't seem to get the right one.

How should I go about converting this code (or any arbitrary code of a similar type) into a custom function?

Answer 1

To construct a generic function that takes a set of equations, variables, limits and plots the output per your code above, you'll need to ensure that the variables (v, ω and t) are properly localized and that global values do not enter your implementation. You can do this with a combination of Block and attributes (ignore the red):

ClearAll@TrajectoryMap
SetAttributes[TrajectoryMap, HoldAll]
TrajectoryMap[
    {eqns__}, 
    {{xvar_, xmin_, xmax_}, {yvar_, ymin_, ymax_}}, 
    {{t_, init_}, tmin_, tmax_}] := 

Block[{t, xvar, yvar, a, b, Pd, PVec},
    Pd = ParametricNDSolve[{eqns, xvar[init] == a, yvar[init] == b}, 
        {xvar[t], yvar[t]}, {t, tmin, tmax}, {a, b}];
    PVec[x_, y_] := {D[(xvar[t] /. Pd)[x, y], t] /. t -> init, 
        D[(yvar[t] /. Pd)[x, y], t] /. t -> init};
    VectorPlot[PVec[x, y], {x, xmin, xmax}, {y, ymin, ymax}]
]

Now you can use this for any set of equations and variables as (using your example):

TrajectoryMap[
    {v'[t] == (v[t] - v[t]^3/3 - ω[t])/0.2, ω'[t] == (1 + v[t])/0.2}, 
    {{v, -2, 1}, {ω, -1.2, 0.5}}, {{t, 0}, 0, 20}]

which is what your example produces.