We present a framework for recovering/approximating unknown time-dependent partial differential equation (PDE) using its solution data. Instead of identifying the terms in the underlying PDE, we seek to approximate the evolution operator of the underlying PDE numerically. The evolution operator of the PDE, defined in infinite-dimensional space, maps the solution from a current time to a future time and completely characterizes the solution evolution of the underlying unknown PDE. Our recovery strategy relies on approximation of the evolution operator in a properly defined modal space, i.e., generalized Fourier space, in order to reduce the problem to finite dimensions. The finite dimensional approximation is then accomplished by training a deep neural network structure, which is based on residual network (ResNet), using the given data. Error analysis is provided to illustrate the predictive accuracy of the proposed method. A set of examples of different types of PDEs, including inviscid Burgers' equation that develops discontinuity in its solution, are presented to demonstrate the effectiveness of the proposed method.ing studies are relatively limited, as they mostly focus on learning certain types of PDE or identifying the exact terms in the PDE from a (large) dictionary of possible terms. The specific novelty of this paper is that the proposed method seeks to recover/approximate the evolution operator of the underlying unknown PDE and is applicable for general class of PDEs. The evolution operator completely characterizes the time evolution of the solution. Its recovery allows one to conduct prediction of the underlying PDE and is effectively equivalent to the recovery of the equation. This is an extension of the equation recovery work from [17], where the flow map of the underlying unknown dynamical system is the goal of recovery. Unlike the ODE systems considered in [17], PDE systems, which is the focus of this paper, are of infinite dimension. In order to cope with infinite dimension, our method first reduces the problem into finite dimensions by utilizing a properly chosen modal space, i.e., generalized Fourier space. The equation recovery task is then transformed into recovery of the generalized Fourier coefficients, which follow a finite dimensional dynamical system. The approximation of the finite dimensional evolution operator of the reduced system is then carried out by using deep neural network, particularly the residual network (ResNet) which has been shown to be particularly suitable for this task [17]. One of the advantages of the proposed method is that, by focusing on evolution operator, it eliminates the need for time derivatives data of the state variables. Time derivative data, often required by many existing methods, are difficult to acquire in practice and susceptible to (additional) errors when computed numerically. Moreover, the proposed method can cope with solution data that are more sparsely or unevenly distributed in time. Since the proposed framework is rather general, we present seve...