What's wrong with SetDelayed? The body of the function works alone, while setdelayed to it as $f[x\_]:=func$ gives the wrong answer

Question

So I was writing a function to get the poincare map of a given dynamic system, for example, the duffing oscillator.

I obtained the time series in advance, then take the poincare map based on this time series.

(All the code is based on this article)

(*define the vector field*)
Duffing[f_, w_] := {v, -0.5*v - x^3 + x + f*Cos[z], w};

(*Obtain the orbits (numerical iteration)*)
OrbitFlow[flow_, vars_, x0_, {t0_, t1_}, opts___] := 
 Module[{rules, eq, l = Length[flow], tf, rhs}, 
  rules = Table[vars[[i]] -> vars[[i]][t], {i, 1, l}];
  rhs = flow /. rules;
  eq = Join[Table[D[vars[[i]] /. rules, t] == rhs[[i]], {i, 1, l}], 
    Table[vars[[i]][0] == x0[[i]], {i, 1, l}]];
  tf = Table[vars[[i]][t], {i, 1, l}];
  sol = NDSolve[eq, tf, {t, t0, t1}, opts]]
(*solving the equation*)

T=2 \[Pi]
sol = OrbitFlow[Duffing[0.39, 1], {x, v, z}, {0, 0, 0}, {0, 1000*T}];

(*to get the set by poincare map*)
poincare[timeseries_,w_]:=Module[{T=2\[Pi]/w,xsol,vsol},
xsol[t_]=x[t]/.timeseries;
vsol[t_]=v[t]/.timeseries;
((Join[xsol[#],vsol[#]]//Evaluate)&) /@Range[800*T,1000*T,T]
];

it doesn't evaluate, and returns:

In:= poincare[sol, 1] 
Out:= {{x[1600 \[Pi]], v[1600 \[Pi]]},...}

instead, if I run the code without function-defined 'poincare', it gives evaluated answers,

In:= Module[{T = 2 \[Pi]/1, xsol, vsol}, xsol[t_] = x[t] /. sol;
 vsol[t_] = v[t] /. sol;
 ((Join[xsol[#], vsol[#]] // Evaluate) &) /@ Range[800*T, 1000*T, T]]
Out:= {{-0.508376, 0.294264},...}

Though I know there are plenty of ways to do the same thing, I just wondered why this is wrong?

---

so it might be something wrong with my settings in NDSolve, in my previous code, I used

NDSolve[eq,{x[t],v[t],z[t]},{~}]

so the output is a list of rules, like:

{{x[t] -> InterpolatingFunction[~],
  v[t] -> InterpolatingFunction[~],
  z[t] -> InterpolatingFunction[~],
}}

so that evaluation at a specific point may not work, although set another function may solve this but yet cause other problems,

In:= x[0.2]/.sol
Out:= {x[0.2]}

In:= temp[t_] = x[t] /. sol; temp[0.2]
Out:= {0.00754592}

Change the variables in NDSolve may solve this, as

NDSolve[eq,{x,v,z},{~}]

then the solution can be evaluated at some point directly

In:= x[0.2]/.sol
Out:= {0.00754592}

But still, I don't know why this could happen that when setdelayed to it, it will not work.

---

Answer 1

This is actually nontrivial: it has to do with how Mathematica manages scope after substituting in values that match a pattern (as in applying :=).

By inspecting poincare[sol, 1] // TracePrint, we see that after substitution, Mathematica actually uses an expression of the following form:

Module[{T$=2 \[Pi],xsol$,vsol$},
 xsol$[t$_] = x[t$] /. {{x[t]->InterpolatingFunction[<<*>>][t], v[t]->InterpolatingFunction[<<*>>][t], z[t]->InterpolatingFunction[<<*>>][t]}};
 vsol$[t$_] = v[t$] /. {{x[t]->InterpolatingFunction[<<*>>][t], v[t]->InterpolatingFunction[<<*>>][t], z[t]->InterpolatingFunction[<<*>>][t]}};
(Evaluate[Join[xsol$[#1],vsol$[#1]]] &) /@ Range[800 T$,1000 T$,T$]]

Mathematica has tried to append $ to the symbol names in an attempt to avoid the very coincidence that you were counting on!

The easiest solution for your case is to extract the InterpolatingFunctions directly as the Head of InterpolatingFunction[<<*>>][t], and set those to xsol and vsol after composing with List (to recover the singleton list structure).

(I'm also not sure why you have Evaluate there, so I took it out, and it seems to work.)

poincare[timeseries_, w_] := 
 Module[{T = 2 \[Pi]/w, xsol, vsol}, 
  xsol = List@*Head[First[x[t] /. timeseries]]; 
  vsol = List@*Head[First[v[t] /. timeseries]];
  (Join[xsol[#], vsol[#]]) & /@ Range[800*T, 1000*T, T]]

(There are other more in-depth ways to avoid those scoping conflicts, but given that this is not essential to what you're doing, I think it makes more sense to just go directly for what you want! But let me know if you're interested!)