Hybrid energy harvesting systems are broadly applied in various fields due to the advantage of improving energy harvesting efficiency. In actual environment, there are many complex phenomena exhibiting jump, flights, rare transition features, and intermittent features, which can be described by systems subjected to non-Gaussian Levy process. Sometimes, it is difficult to build mathematical models of complex hybrid energy harvesting systems precisely, especially for those driven by non-Gaussian Levy noise. With the development of simulation capabilities and observing techniques recently, massive noise measurement data or simulating data can be feasibly obtained and there are many existing techniques devoted to discovering governing laws from abundant data. In this paper, we aim to extract the system equations from observed data influenced by non-Gaussian Levy noise via using a data-driven method. The expressions of drift term, diffusion term and Levy term can be approximated with the help of Fokker-Planck equation and non-local Kramers-Moyal formulae, and the coefficients of the expressions are learned by utilizing a sparse regression approach in the least square sense. Three examples are given to demonstrate the effectiveness of the method. Results show that the approach can well be applied to not only hybrid energy harvesting systems under Gaussian Brownian process but also the systems subjected to non-Gaussian Levy process. Additionally, the relations between the demarcation parameter and an indicator denoted as Ratio for different time steps are analyzed, and results demonstrate that the indicator can be regarded as the criterion of selecting the appropriate demarcation parameter.