We address the problem of inferring an undirected graph from nodal observations, which are modeled as non-stationary graph signals generated by local diffusion dynamics that depend on the structure of the unknown network. Using the so-called graph-shift operator (GSO), which is a matrix representation of the graph, we first identify the eigenvectors of the shift matrix from observations of the diffused signals, and then estimate the eigenvalues by imposing desirable properties on the graph to be recovered. Different from the stationary setting where the eigenvectors can be obtained directly from the covariance matrix of the measurements, here we need to estimate first the unknown diffusion (graph) filter - a polynomial in the GSO that preserves the sought eigenbasis. To carry out this initial system identification step, we exploit different sources of information on the arbitrarily-correlated input signal driving the diffusion on the graph. We first explore the setting where the observations, the input information, and the unknown graph filter are linearly related. We then address the case where the relation is given by a system of matrix quadratic equations, which arises in pragmatic scenarios where only the second-order statistics of the inputs are available. While such a quadratic filter identification problem boils down to a non-convex fourth-order polynomial minimization, we discuss identifiability conditions, propose algorithms to approximate the solution, and analyze their performance. Numerical tests illustrate the effectiveness of the proposed topology inference algorithms in recovering brain, social, financial, and urban transportation networks using synthetic and real-world signals.