Finite-width fully connected neural networks with Gaussian initialization deviate from their infinite-width Gaussian limit through non-vanishing higher-order cumulants. In this talk, I present multidimensional Edgeworth expansions of arbitrary order for neural network outputs evaluated on a finite collection of inputs, providing a systematic way to approximate these non-Gaussian effects. Under the assumption that the limiting Gaussian covariance matrix is invertible and that the activation function is polynomially bounded, we obtain upper bounds of order $n^{-m}$ in total variation distance between the true network law and its Edgeworth approximation of order $4m-2$, together with matching lower bounds. Beyond neural networks, the results apply to general sequences of conditionally Gaussian vectors converging to a non-degenerate Gaussian limit. As an application, I discuss quantitative bounds for Bayesian neural networks, measuring the error introduced when replacing the prior distribution with its Edgeworth approximation.