Scientific machine learning has been successfully applied to inverse problems and PDE discoveries in computational physics. One caveat of current methods however is the need for large amounts of (clean) data in order to recover full system responses or underlying physical models. Bayesian methods may be particularly promising to overcome these challenges as they are naturally less sensitive to sparse and noisy data. In this paper, we propose to use Bayesian neural networks (BNN) in order to: 1) Recover the full system states from measurement data (e.g. temperature, velocity field, etc.). We use Hamiltonian Monte-Carlo to sample the posterior distribution of a deep and dense BNN, and show that it is possible to accurately capture physics of varying complexity without overfitting. 2) Recover the parameters in the underlying partial differential equation (PDE) governing the physical system. Using the trained BNN as a surrogate of the system response, we generate datasets of derivatives potentially comprising the latent PDE of the observed system and perform a Bayesian linear regression (BLR) between the successive derivatives in space and time to recover the original PDE parameters. We take advantage of the confidence intervals on the BNN outputs and introduce the spatial derivative variance into the BLR likelihood to discard the influence of highly uncertain surrogate data points, which allows for more accurate parameter discovery. We demonstrate our approach on a handful of example applied to physics and non-linear dynamics.