# Writing NDSolve-like functions

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.

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?

-
If I understood it correctly, you want TrajectoryMap[...] to do what is being done in the code block above, right? If so, where do you specify the ICs and the plot ranges? –  rm -rf Oct 29 '13 at 20:55
@rm-rf The ICs are the ParametricNDSolve parameters, they get substituted with x,y later on. The vector plot is a plot of $(v'(t),\omega'(t))$ at various initial conditions. The plot ranges are specified in the {{v[t],-2,1}, {ω[t],-1.2,0.5}} (i.e. the {{var1,varmin,varmax},{var2,varmin,varmax}}) part. –  Manishearth Oct 29 '13 at 20:59
Yes, but for a generic function, shouldn't you allow for an IC point other than 0 and t other than 0–20? –  rm -rf Oct 29 '13 at 21:04
@rm-rf Well, if that can be genericized I don't mind, but when I use this set of code I usually don't need much of a range for t -- in fact, none at all. I'll add that bit, though. –  Manishearth Oct 29 '13 at 21:06

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.

-
Thanks! I didn't know about SetAttributes, (or Block, I've been using Module which is lexical.), which must be where I went wrong. I also didn't know about the ability to use lists in the function declaration itself. I think that there may be a syntax error here, though, the t,xvar,yvar shouldn't be in the Block statement. Not sure. –  Manishearth Oct 30 '13 at 2:28
Nope, not a syntax error. As I said in my answer, ignore the red :) If you want to see why, remove them from Block and redefine TrajectoryMap, evaluate t = v = ω = 1;, followed by the TrajectoryMap[...] line. –  rm -rf Oct 30 '13 at 2:37
It seems to work both ways, though. What's the difference? Is it preventing the variables from messing with outside variables? –  Manishearth Oct 30 '13 at 2:40
@Manishearth I updated my comment while you were writing yours... please see my previous comment. Yes, it is what you said. –  rm -rf Oct 30 '13 at 2:41
Ah, thanks. Why does it work in this case though? When I try t = 5; Print[t]; Block[{t}, Print[t]; t = 6;]; Print[t], the t in the block is taken as a fresh variable. Whereas here it seems to be substituted in from the outside. –  Manishearth Oct 30 '13 at 2:44