Green's functions are fundamental mathematical objects that serve as the impulse response of a linear differential operator, providing deep mechanistic insight into the physical systems they describe. This work presents a data-driven framework to learn and interpolate these vital solution operators directly from observational data, where the governing Partial Differential Equations of the system are unknown. The first work introduces two methods for learning a discrete, low-rank approximation of a system's Green's function from pairs of system excitations and responses: The first approach uses the proper orthogonal decomposition (POD) modes of the system's output as a surrogate for the Green's function's eigenvectors and subsequently fits the corresponding eigenvalues from data. The second, more accurate method employs a generalization of the randomized singular value decomposition (SVD) to construct the approximation. To handle systems with varying parameters, a principled manifold interpolation scheme was developed. This scheme maps the learned discrete eigenmodes to a tangent space for interpolation and then retracts them back to the manifold, allowing for the prediction of the Green's function at unseen parameter instances. Building upon the initial work, a more advanced, mesh-independent framework was developed to represent the learned Green's function as a continuous, bivariate function. In this approach, a Rational Neural Network first learns the Green's function from data. Subsequently, a Singular Value Expansion (SVE) of this continuous function is constructed using a representation in a Chebyshev basis. This method improves data efficiency and extends the scope of applicability from self-adjoint to non-self-adjoint operators. The manifold interpolation algorithm was then generalized to operate on this continuous representation, enabling the interpolation of Quasimatrices (matrices with functional columns) on an infinite-dimensional analogue of the Stiefel manifold. Numerical results demonstrate that this comprehensive framework learns high-fidelity approximations of solution operators and is robust to significant noise in the training data.