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?

  • $\begingroup$ 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? $\endgroup$
    – rm -rf
    Oct 29, 2013 at 20:55
  • $\begingroup$ @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. $\endgroup$ Oct 29, 2013 at 20:59
  • $\begingroup$ Yes, but for a generic function, shouldn't you allow for an IC point other than 0 and t other than 0–20? $\endgroup$
    – rm -rf
    Oct 29, 2013 at 21:04
  • $\begingroup$ @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. $\endgroup$ Oct 29, 2013 at 21:06

1 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):

SetAttributes[TrajectoryMap, HoldAll]
    {{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):

    {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.

  • $\begingroup$ 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. $\endgroup$ Oct 30, 2013 at 2:28
  • $\begingroup$ 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. $\endgroup$
    – rm -rf
    Oct 30, 2013 at 2:37
  • $\begingroup$ It seems to work both ways, though. What's the difference? Is it preventing the variables from messing with outside variables? $\endgroup$ Oct 30, 2013 at 2:40
  • $\begingroup$ @Manishearth I updated my comment while you were writing yours... please see my previous comment. Yes, it is what you said. $\endgroup$
    – rm -rf
    Oct 30, 2013 at 2:41
  • $\begingroup$ 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. $\endgroup$ Oct 30, 2013 at 2:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.