There's a minor mistake in your question, $\langle T j^\mu(x)\psi(y)\overline{\psi}(z)\rangle$ is not the same diagram as $\langle T A^\mu(x)\psi(y)\overline{\psi}(z)\rangle$, the latter contains one extra bare photon propagator, this can be seen by considering possible contractions of the fields. So your question should be modified as: how to prove the equation
$\langle T A^\mu(x)\psi(y)\overline{\psi}(z)\rangle=\int d^4s G^\mu_{\ \ \nu}(x,s)\langle T j^\nu(s)\psi(y)\overline{\psi}(z)\rangle\cdots\cdots(1),$
where $G(x,s)$ is the bare photon propagator?
Like you said, this can be justified just by inspecting the corresponding Feynman diagrams, thus you already have a perturbative proof! So I guess what you really want is a non-perturbative proof, here it goes:
Recall the bare photon propagator is really just the inverse operator(alias Green's function) of the differential operator $d^{\mu}_{\ \ \nu} $in the field equation acting on $A^\nu$, for example, with a Lorentz gauged Lagrangian $\mathcal{L}=-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}+\frac{1}{2\xi}(\partial_\mu A^\mu)^2-j^\mu A_\mu+\mathcal{L}_{matter}$, the field equation is
$d^{\mu}_{\ \ \nu}A^\nu=[-\delta^\mu_{\ \ \nu}\partial^2+(1+\frac{1}{\xi})\partial^\mu\partial_\nu]A^\nu=j^\mu\cdots\cdots(2),$
and the propagator is just $G^\mu_{\ \ \nu}=(d^{-1})^\mu_{\ \ \nu}$, that is,
$[-\delta^\mu_{\ \ \nu}\partial^2+(1+\frac{1}{\xi})\partial^\mu\partial_\nu]G^\nu_{\ \ \rho}(x,s)=\delta^\mu_{\ \ \rho}\delta^4(x-s)$. So if we "de-convolute"(i.e. act $d^{\mu}_{\ \ \nu}$ on) both sides of equation (1), we would get an equivalent expression as
$d^{\mu}_{x\ \nu}\langle T A^\nu(x)\psi(y)\overline{\psi}(z)\rangle=\langle T j^\mu(x)\psi(y)\overline{\psi}(z)\rangle\cdots\cdots(3),$
where the subscript $x$ in $d^{\mu}_{x\ \nu}$ indicates the space-time variable it's acting on. So we might as well just prove (3). Now (3) looks like (2) a lot, but one has to be careful because $d^{\mu}_{x\ \nu}$ involves time derivatives hence does not commute with the time-ordering symbol. One can expand time-ordering in terms of step functions and differentiate term by term carefully, then summing up together the resulting terms will show (3) is indeed true, as if $d^{\mu}_{x\ \nu}$ commutes with time ordering. This is a lot of work, but one can use a much more efficient method called
Dyson-Schwinger(D-S) equation(although there's still something subtle about it that I don't understand, but it has never disappointed me in practice). Applying D-S equation to $\langle T\psi(y)\overline{\psi}(z)\rangle$ by considering variations in gauge fields will immediately lead to equation (3).
Last but not least(maybe OP is aware of the following, but I can't tell from the question, so I'll write it anyway), I must emphasize in general a simple relation like equation (1) doesn't hold if you have more than one gauge field in the time-ordered product. For example, you can't get $\langle T A^\mu(x) A^\rho(y) j^\sigma(z) \rangle$ just by attaching two bare photon propagators to $\langle T j^\mu(x) j^\rho(y) j^\sigma(z) \rangle$. Again you can convince yourself in two different ways:
In terms of Feynman diagrams, now you can have a Wick contraction between external points $A^\mu(x)$ and $A^\rho(y)$(unlike when you only have one gauge field operator), this is just a bare photon propagator. The differences can be visually captured by the following figures(note that a blob is not necessarily a connected part):
The second way is again in terms of D-S equation, you start from $\langle j^\sigma(z)\rangle$ and apply D-S equation iteratively twice, you'll get exactly the right terms corresponding to the third diagram above, and as it turns out the contraction between $A^\mu(x)$ and $A^\rho(y)$ corresponds to what's called a "contact term" in D-S equation. This time you can't pretend $d^{\mu}_{\ \ \nu}$ commutes with time-ordering symbol!
