Creating a vibrato in Mathematica

How about solving the harmonic oscillator equation with a time-varying frequency?

w[t_] := 2 Pi (440 + 5 Sin[10 * 2 Pi t])

func = NDSolveValue[{y''[t] + w[t]^2 y[t] == 0, y[0] == 1, y'[0] == 0}, y, {t, 0, 4}];

Play[func[t], {t, 0, 4}]

As an alternative, let's let the derivative of the phase vary with time.

A fixed tone will have $d \phi / d t = 2 \pi \cdot 440$. A vibrato that you want should have $d \phi / d t = 2 \pi (440 + \sin (2\pi\cdot5\cdot t))$. Integrating:

Integrate[440*2 π + 2 π Sin[5*2 π x], {x, 0, t}]
(* 880 π t + 2/5 Sin[5 π t]^2 *)

Playing:

Play[Sin[880 π t + 2/5 Sin[5 π t]^2], {t, 0, 4}]