I have a code, which is quite large. There, I consider parameters that are used inside series of procedures and numerically solved differential equations. At the end I plot the solution for those parameters. If I want to consider another situation, I have to consider another set of parameters, and plot them.

I would like to have a plot that includes, say, three cases of the parameters, and plot them together to see how the solutions changes with the paremeter.

My program is divided in blocks, each of them shows some particular results I want in the intermediate steps. So, it is like turning the code into a Table. You could translate my question into putting a DO or some repetitive command of a programming language, so the last result looks like running the code starting from a different parameters each time.

Edited part: the code


In[1]:= \[Epsilon] := 10

In[2]:= L := 10

In[3]:= d := 1

In[4]:= \[Phi] := \[Pi]/4

In[5]:= \[Chi]m := \[Pi]/4

In[6]:= (*\[Chi]p:=1/1000*)

In[7]:= \[Alpha]p := 1

In[8]:= N5 := 1

In[9]:= N3 := 100

In[10]:= \[CapitalDelta]N3 := 2

In[11]:= gYM := 1

In[12]:= \[Delta] = 
  1/2 Log[1/(
     gYM^2 N5^2 (2 N3 - \[CapitalDelta]N3)) (2 gYM^2 N3 N5^2 + 
       4 \[Pi]^2 \[CapitalDelta]N3^2 + 
       Sqrt[(2 gYM^2 N3 N5^2 + 4 \[Pi]^2 \[CapitalDelta]N3^2)^2 - 
         gYM^4 N5^4 (4 N3^2 - \[CapitalDelta]N3^2)])];

In[13]:= \[Alpha] = -(N5/4) Cosh[\[Delta]] + 
   Sqrt[(\[Pi]^2 N3)/gYM^2 + N5^2/16 Cosh[\[Delta]]^2];

In[14]:= \[Alpha]h = (gYM^2 \[Alpha])/(4 \[Pi]);

(*harmonic functions and definitions*)

In[15]:= h1 = \[Alpha]p (-I \[Alpha] Sinh[v] - 
      N5/4 Log[Tanh[(I \[Pi])/4 - (v - \[Delta])/
          2]]) + \[Alpha]p (I \[Alpha] Sinh[vb] - 
      N5/4 Log[Tanh[-((I \[Pi])/4) - (vb - \[Delta])/2]]);

In[16]:= h2 = \[Alpha]p \[Alpha]h (Cosh[v] + Cosh[vb]);

In[17]:= w = D[D[h1 h2, vb], v];

In[18]:= F1 = 2 h1 h2 D[h1, v] D[h1, vb] - h1^2 w;

In[19]:= F2 = 2 h1 h2 D[h2, v] D[h2, vb] - h2^2 w;

In[20]:= f42s = 2 ((F1 F2)/w^2)^(1/4);

In[21]:= \[Rho]2s = 2 h2^2 ((F1 w^2)/F2^3)^(1/4);

In[22]:= subv = {v -> x[x2] + I y[x2], vb -> x[x2] - I y[x2]};

Coefficient functions inside the eoms


In[24]:= (*these functions will be evaluated at the initial values given below, BUT \
they give undefined values or non-significant digits depending on the \

In[25]:= Logfx = D[Log[f42s /. subv], x[x2]];

In[26]:= Logfy = D[Log[f42s /. subv], y[x2]];

In[27]:= Log\[Rho]x = D[Log[\[Rho]2s /. subv], x[x2]];

In[28]:= Log\[Rho]y = D[Log[\[Rho]2s /. subv], y[x2]];

In[29]:= f\[Rho]x = D[f42s /. subv, x[x2]]/(\[Rho]2s /. subv);

In[30]:= f\[Rho]y = D[f42s /. subv, y[x2]]/(\[Rho]2s /. subv);

(Lagrangian and equations of motion*)

In[31]:= (*f42 is Subscript[f, 4](x,y)^2 and \[Rho]2 is \[Rho](x,y)^2*)

In[32]:= (*Subscript[f, 4](x,y) and \[Rho](x,y) are given*)


In[34]:= subeq = { f42[x[x2], y[x2]] -> f42, 
RowBox[{"0", ",", "1"}], ")"}],
MultilineFunction->None]\)[x[x2], y[x2]] -> f42y, 
RowBox[{"1", ",", "0"}], ")"}],
MultilineFunction->None]\)[x[x2], y[x2]] -> 
    f42x, \[Rho]2[x[x2], y[x2]] -> \[Rho]2, 
RowBox[{"1", ",", "0"}], ")"}],
MultilineFunction->None]\)[x[x2], y[x2]] -> \[Rho]2x, 
RowBox[{"0", ",", "1"}], ")"}],
MultilineFunction->None]\)[x[x2], y[x2]] -> \[Rho]2y};

In[35]:= subu = {f42x -> A f42, f42y -> B f42};

In[36]:= subxy = {f42x -> F \[Rho]2, 
   f42y -> G \[Rho]2, \[Rho]2x -> H \[Rho]2, \[Rho]2y -> J \[Rho]2};

In[37]:= subwarp2 = {A -> Logfx, B -> Logfy, H -> Log\[Rho]x, J -> Log\[Rho]y, 
   F -> f\[Rho]x, G -> f\[Rho]y};

In[38]:= subwarphf = {A -> HoldForm@Logfx, B -> HoldForm@Logfy, 
  H -> HoldForm@Log\[Rho]x, J -> HoldForm@Log\[Rho]y, F -> HoldForm@f\[Rho]x, 
  G -> HoldForm@f\[Rho]y}

Out[38]= {A -> \!\(\*
HoldForm]\), B -> \!\(\*
HoldForm]\), H -> \!\(\*
HoldForm]\), J -> \!\(\*
HoldForm]\), F -> \!\(\*
HoldForm]\), G -> \!\(\*

In[39]:= (*Lagrangian*)

In[40]:= Lag = f42[x[x2], y[x2]] (u'[x2]^2/u[x2]^2 + 2/u[x2]^2) + \[Rho]2[x[x2], 
     y[x2]] (x'[x2]^2 + y'[x2]^2);

In[41]:= (*SET OF EQUATIONS OF MOTION where A,B,F,G,H and J are defined above*)

In[42]:= (*eom for u(x2)*)

In[43]:= equ =  (u[x2]^2)/(2 f42) (D[Lag, u[x2]] - D[D[Lag, u'[x2]], x2]) /. subeq /. 
    subu // Expand;

In[44]:= equ /. subwarphf

Out[44]= -(2/u[x2]) + Derivative[1][u][x2]^2/u[x2] - \!\(\*
HoldForm]\) Derivative[1][u][x2] Derivative[1][x][x2] - \!\(\*
HoldForm]\) Derivative[1][u][x2] Derivative[1][y][x2] - (u^\[Prime]\[Prime])[

In[45]:= (*eom for x(x2)*)

In[46]:= eqx =  1/(2 \[Rho]2) (D[Lag, x[x2]] - D[D[Lag, x'[x2]], x2]) /. subeq /. 
    subxy // Expand;

In[47]:= eqx /. subwarphf

Out[47]= \!\(
HoldForm]\)/u[x2]^2 + (\!\(\*
HoldForm]\) Derivative[1][u][x2]^2)/(2 u[x2]^2) - 1/2 \!\(\*
HoldForm]\) Derivative[1][x][x2]^2 - \!\(\*
HoldForm]\) Derivative[1][x][x2] Derivative[1][y][x2] + 1/2 \!\(\*
HoldForm]\) Derivative[1][y][x2]^2 - (x^\[Prime]\[Prime])[x2]

In[48]:= (*eom for y(x2)*)

In[49]:= eqy =  1/(2 \[Rho]2) (D[Lag, y[x2]] - D[D[Lag, y'[x2]], x2]) /. subeq /. 
    subxy // Expand;

In[50]:= eqy /. subwarphf

Out[50]= \!\(
HoldForm]\)/u[x2]^2 + (\!\(\*
HoldForm]\) Derivative[1][u][x2]^2)/(2 u[x2]^2) + 1/2 \!\(\*
HoldForm]\) Derivative[1][x][x2]^2 - \!\(\*
HoldForm]\) Derivative[1][x][x2] Derivative[1][y][x2] - 1/2 \!\(\*
HoldForm]\) Derivative[1][y][x2]^2 - (y^\[Prime]\[Prime])[x2]

In[51]:= (*collecting the equations*)

In[52]:= pdes = {equ == 0, eqx == 0, eqy == 0} /. subwarp2;

(*Boundary conditions*)

In[53]:= (*initial and final values of x2*)

In[54]:= x20 = -d Cos[\[Phi]];

In[55]:= x21 = d Cos[\[Phi]];

In[56]:= (*initial values of u,x and y*)

In[57]:= u0 = \[Epsilon] Sqrt[1 + ((L - d Sin[\[Phi]])/\[Epsilon])^2];

In[58]:= x0 = ArcSinh[(L - d Sin[\[Phi]])/\[Epsilon]];

In[59]:= y0 = \[Pi]/2 - \[Chi]m;

In[60]:= (*Final values*)

In[61]:= u1 = \[Epsilon] Sqrt[1 + ((L + d Sin[\[Phi]])/\[Epsilon])^2];

In[62]:= x1 = ArcSinh[(L + d Sin[\[Phi]])/\[Epsilon]];

In[63]:= y1 = \[Pi]/2 - \[Chi]p;

In[64]:= (*required condition*)

In[65]:= Ld = L + d Sin[\[Phi]] // N

Out[65]= 10.7071

In[66]:= L - d Sin[\[Phi]] > 0

Out[66]= True

In[67]:= (*boundary conditions*)

In[68]:= in = {u0, x0, y0} // N

Out[68]= {13.6513, 0.83048, 0.785398}

In[69]:= out = {u1, x1, y1} // N

Out[69]= {14.6507, 0.9305, 1.5708 - 1. \[Chi]p}

In[70]:= (*bvp*)

In[71]:= {x20, x21} // N

Out[71]= {-0.707107, 0.707107}

In[72]:= x21 - x20 // N

Out[72]= 1.41421

In[73]:= slu = (u1 - u0)/(x21 - x20) // N

Out[73]= 0.706665

In[74]:= slx = (x1 - x0)/(x21 - x20) // N

Out[74]= 0.0707254

In[75]:= sly = (y1 - y0)/(x21 - x20) // N

Out[75]= 0.707107 (0.785398 - 1. \[Chi]p)

In[76]:= bcs = {x[x20] == x0, u[x20] == u0, y[x20] == y0, x[x21] == x1, u[x21] == u1, 
   y[x21] == y1};

In[77]:= (*ibv*)

In[78]:= bcsin = {u[x20] == u0, x[x20] == x0, y[x20] == y0};

In[79]:= (*values+parameter*)

In[80]:= (*bcsp={bcsin,u'[x20]\[Equal]0.3504+a,x'[x20]\[Equal]-18.9337+b,y'[x20]\

In[81]:= (*first derivatives*)

In[82]:= bcspg = {bcsin, u'[x20] == r, x'[x20] == s, y'[x20] == t} // Flatten;

(*Solving eoms*)

In[83]:= solspg = ParametricNDSolveValue[{pdes, bcspg}, {u[x2], x[x2], y[x2]}, {x2, 
    x20, x21}, {r, s, t}];

In[84]:= (*Plot[Evaluate[Table[solspg[r,s,t],{r,0,1,0.1},{s,-20,10,1},{t,15,25,1}]],{\

In[85]:= out

Out[85]= {14.6507, 0.9305, 1.5708 - 1. \[Chi]p}

In[86]:= outn = solspg[slu, slx, 0] /. x2 -> x21 // Chop

Out[86]= {14.4611, 0.963592, 0.777848}

In[87]:= Solve[{(outn[[1]])^2 == \[Epsilon]n^2 + Ldn^2, 
   outn[[2]] == ArcSinh[Ldn/\[Epsilon]n]}, {\[Epsilon]n, Ldn}] // Last

Out[87]= {\[Epsilon]n -> 9.63237, Ldn -> 10.7862}

In[88]:= Solve[outn[[3]] == y1, \[Chi]p]

Out[88]= {{\[Chi]p -> 0.792948}}

In[89]:= \[Chi]m // N

Out[89]= 0.785398

In[90]:= Ld

Out[90]= 10.7071

In[91]:= Plot[Evaluate[solspg[slu, slx, 0][[1]]], {x2, x20, x21}, ImageSize -> Large, 
 LabelStyle -> Directive[Medium], AxesLabel -> {"x2", "u"}, 
 PlotRange -> {{x20, x21}, All}]

In[92]:= Plot[Evaluate[solspg[slu, slx, 0][[2]]], {x2, x20, x21}, ImageSize -> Large, 
 LabelStyle -> Directive[Medium], AxesLabel -> {"x2", "u"}, 
 PlotRange -> {{x20, x21}, All}]

In[93]:= Plot[Evaluate[solspg[slu, slx, 0][[3]]], {x2, x20, x21}, ImageSize -> Large, 
 LabelStyle -> Directive[Medium], AxesLabel -> {"x2", "u"}, 
 PlotRange -> {{x20, x21}, All}]

Here you see a set of parameters, some of them I want to fix and some others I want to change and compare solutions in one single plot. As mentioned, my code is given in blocks to show some steps in the middle, some useful numbers, etc. Please, forget the (*comments*) lines. The parameters I want to manipulate are N5, N3 and \[CapitalDelta]N3. In the program below I consider N5=1, N3=100 and \[CapitalDelta]N3=2, but I'd like to modify these valuesand compare the cases N5=1,10,100, N3=10,100 and \[CapitalDelta]N3=1,10,100, and see if the differential equations have solution, etc. Everything without changing the structure of the code, the blocks.

  • $\begingroup$ If you want to execute a chunk of codes repeatedly, I think the best way is to make them a function. $\endgroup$ Commented Dec 6, 2017 at 5:00
  • $\begingroup$ @AlexanderZeng Yes, but my question is about of considering some line or command at the beginning of the code only that allows to evaluate several cases in one single run. This is because my code is given inn blocks, etc. But yes, that would be an option. $\endgroup$ Commented Dec 6, 2017 at 13:57
  • $\begingroup$ It might be better if you provide a MWE. $\endgroup$ Commented Dec 6, 2017 at 14:05
  • $\begingroup$ Without more detailed information, I can only describe a line of thought. First, just in my opinion, it is better to make the two things independent: obtaining data results and making figures. So at the beginning you have several groups of parameters para = {p1, p2, ..., pn} (each pi can be a List.), for each group you want to get the corresponding data for plotting dataTotal = {data1, data2, ..., datan}. So a big function f can be defined to have the mapping datai=f[pi]. Then you can have dataTotal = f@@@para, after which you can use it for plotting. $\endgroup$ Commented Dec 6, 2017 at 14:16
  • $\begingroup$ @AlexanderZeng Code added. Thanks for your help. $\endgroup$ Commented Dec 6, 2017 at 17:22

2 Answers 2


Module can be used to turn sequential code into a function. As a simple example, suppose I use the following code to make a plot of $A\sin(kx + \delta)$, where (for whatever reason) I want $A = kd$:

k = 1;
delta = 2;
A = k*delta;
Plot[ A * Sin[ k*x + delta], {x, 0, 4*Pi}]

I can use Module to define a function that takes values of k and delta and spits out a plot of $\sin(kx + \delta)$:

sinplot[kin_, deltain_] := Module[{k = kin, delta = deltain, A},
    A = k * delta;
    Plot[ A * Sin[ k*x + delta], {x, 0, 4*Pi}]

Invoking sinplot[1, 2] will then give me the same plot as I had in the previous code block.

As far as putting multiple plots together, that's a job for Show:

plot1 = sinplot[1,2];
plot2 = sinplot[2,2];
plot3 = sinplot[2,1];
Show[plot1, plot2, plot3]

enter image description here

However, the combined plot from Show "inherits" its properties from the individual plots, which means that they'll generally all have the same attributes. A better solution might be to have your Module output the numerical solution to the ODE rather than the plot, and then use something like Plot[{soln1, soln2, soln3}] instead.

  • $\begingroup$ Thanks. I had some problems adding the Module you suggested. I dont know if this is because the block structure of the code, or because there are PDEs in the middle steps that need to be solved with numerical values. I added the code, please give a look if you want. $\endgroup$ Commented Dec 6, 2017 at 17:24

For completeness, be aware that this kind of code works on Mathematica :

A = k*delta;
Plot[ A * Sin[ k*x + delta], {x, 0, 4*Pi}])

k = 1;     delta = 2    ;automatisedCode,
k = 2    (*delta = 2;*) ;automatisedCode,
(*k = 2;*) delta = 1    ;automatisedCode

enter image description here

No need to use Module[...]. Module is intended to localize some variables, it is not a kind of box for automatised-code embedding (function, module, routine etc... in other languages).

Of course, this kind of code is not recommended in large projects.

It's interesting to use automatisedCode[] instead of automatisedCode, specially if you have in mind to pass afterwards k or delta as parameter of automatisedCode[] (by using automatisedCode[k]). If you have used automatisedCode without the [] (even only once), you have to do do Clear[f] before defining automatisedCode[k]. (if you don't do that, it may cause subtle issues)

  • 1
    $\begingroup$ Better to use DownValues than OwnValues, I think. Procedurally the difference is minimal, but delayed OwnValues can have a way of causing subtle issues, so unless it makes more conceptual sense (or I want the OwnValue to be overwrite-able) I tend to stick to DownValues. $\endgroup$
    – b3m2a1
    Commented Dec 5, 2017 at 17:16
  • $\begingroup$ @b3m2a1 M. Wizard had the same reaction in another answer. See my text here : mathematica.stackexchange.com/questions/134505/… $\endgroup$
    – andre314
    Commented Dec 5, 2017 at 18:20
  • $\begingroup$ @b3m2a1 & Co I have added a paragraph at the end of my answer that takes account of your comment. I mean Ownvalue corresponds to use automatisedCode, Downvalue corresponds to use automatisedCode[] $\endgroup$
    – andre314
    Commented Dec 6, 2017 at 2:37
  • $\begingroup$ @andre & Co, suggestions give problems when there are PDEs inside the code. $\endgroup$ Commented Dec 8, 2017 at 0:50

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.