In this work, a novel approach in determining the first- and second-order frequency-domain Volterra kernels for weakly nonlinear partial differential equations (PDEs) with a forcing term in semidiscrete form based on the application of the harmonic probing method is presented. This represents a formal extension of the linearized frequency-domain (LFD) methods to a nonlinear framework, leading to a so-called LFD method for a second-order functional series expansion method. The method allows for the representation of weak nonlinearities by solving two input-independent linear algebraic systems of equations in the frequency domain, and thus circumvents the solution of the nonlinear PDE by numerical integration for each different input, representing a nonlinear reduced-order model (ROM) for the physical phenomena. The general form of the equations is derived, and several applications are discussed. The proposed method overcomes two constraints present in other methods for the solution of nonlinear PDEs, namely, the consideration of exclusively periodic solutions as in the harmonic balance method and the dependency of the kernels with the input signal as in the Volterra kernel identification methods.