I wouldn't be surprised if someone wrote a tool that allows you to use a reference to an equation as the function to plot :).

You can write (or derive) the expression using sympy, then have sympy generate a numpy expression that can be evaluated. Sympy can also generate the LaTeX code for any expression. So while that isn't an in-LaTeX solution, it may be close to what you want.

Johansson's "Numerical Python" shows several examples of this. I will scavenge one of his examples below (trusting it falls under "fair use", and hoping I transcribe it correctly -- note I have left out the imports). The example uses sympy to generate and plot Taylor series expansions of sin(x).

The key bit to look for in the example is `sympy.lambdify()`.

    sym_x = sympy.Symbol("x")
    x = np.linspace(-2 * np.pi, 2 * np.pi, 100)

    def sin_expansion(x, n):
        return sympy.lambdify(sym_x, sympy.sin(sym_x).series(n=n+1).removeO(), 'numpy')(x)

    fig, ax = plt.subplots()
    ax.plot(x, np.sin(x), linewidth=4, color="red", label='exact')
    colors = ["blue", "black"]
    linestyles = [':', '-.', '--']

    for idx, n in enumerate(range(1, 12, 2)):
        ax.plot(x, sin_expansion(x, n), color=colors[idx // 3],
            linestyle=linestyles[idx % 3], linewidth=3,
            label="order %d approx." % (n+1))

    ax.set_ylim(-1.1, 1.1)
    ax.set_xlim(-1.5*np.pi, 1.5*np.pi)

    ax.legend(bbox_to_anchor=(1.02, 1), loc=2, borderaxespad=0.0)
    fig.subplots_adjust(right=.75)

Also, (v)phantom and smash are something I really should have learned before all those fancy packages, nowadays I'm mostly using context anyways.