# How can I change the scoping of a pure function returned by a function in a package?

I have written a small Mathematica package that contains a module that solves some differential equation with initial conditions. Suppose that the solution to this differential equation defines some other functions y(t) and z(t). I want this function to return these functions y and z as rules, much as NDSolve returns a rule for evaluating a function.

Here is my first naive attempt to do this:

SolveEquation[x_, xp0_, x0_] := Module[{y, z},
sol = Flatten[
NDSolve[{x''[t] + x'[t]^3 == x[t]^2, x'[0] == xp0,
x[0] == x0}, {x}, {t, 0, 100}]];
y[t_] := (x[t] /. sol)^1/2;
z[t_] := y'[t] /. y[t];
Return[Join[sol, {Symbol["y"] -> Function[t, y[t]],Symbol["z"] -> Function[t, z[t]]}]]]

I can then correctly access the values of these functions with
rules=SolveEquation[x,0,100];
y[3]/.rules
9.89546

However, a call to rules itself gives me:
{x -> InterpolatingFunction[{{0.,100.}},<>],y -> Function[t$, y$5338[t$]], z -> Function[t$, z$5338[t$]]}

This is okay but in terms of human readability (the code will be used by individuals besides myself), I would much prefer:
{x -> InterpolatingFunction[{{0.,100.}},<>],y -> Function[t, y[t]],
z -> Function[t,z[t]}

or even better:
{x -> InterpolatingFunction[{{0.,100.}},<>],y -> x[t]^1/2,
z -> y'[t]/y[t]}

or some other way in which the dependence of each function on the others is made explicit.

Is it possible to achieve this result by somehow changing the scoping of these variables or do I need to just settle for my method that works but appears sloppy?

-
Note that in NDSolve the variables that are returned in the rule are actually part of the NDSolve function call. In other words something like myF[a_, x_Symbol] := x -> a + 1 –  user21 Jul 11 '14 at 19:30
You are correct. I am asking if it is possible to replicate this outcome without having to accept an argument for each of the 20 functions that I need to define. –  Jerro39 Jul 11 '14 at 19:48
A function with 20 arguments, I think, is very difficult to use. If you really want to go that way, you could use a Format for hide the ugly stuff. –  user21 Jul 11 '14 at 19:55

This could be a way of achieving what you want (if I understood you correctly):

ClearAll[SolveEquation];
SolveEquation[x_Symbol, xp0_, x0_] :=
Block[{x, y, z},
Module[{sol},
sol = First@
Flatten[NDSolveValue[{x''[t] + x'[t]^3 == x[t]^2, x'[0] == xp0,
x[0] == x0}, {x}, {t, 0, 100}]];
{
x -> sol,
y -> (x[#]^(1/2) &),
z -> (y'[#]/y[#]] &)
}
]
]


Use ReplaceRepeated for evaluation:

y'[3] //. rules


And rules returns:

{x -> InterpolatingFunction[{{0.,100.}},"<>"],
y -> (Sqrt[x[#1]] &),
z -> (y'[#1]/y[#1] &)}

-
Thank you! This is exactly what I was looking for. –  Jerro39 Jul 11 '14 at 21:55