In this work, we propose a robust Bayesian sparse learning algorithm based on Bayesian group Lasso with spike and slab priors for the discovery of partial differential equations with variable coefficients. Using the samples draw from the posterior distribution with a Gibbs sampler, we are able to estimate the values of coefficients, together with their standard errors and confidence intervals. Apart from constructing the error bars, uncertainty quantification can also be employed for designing new criteria of model selection and threshold setting. This enables our method more adjustable and robust in learning equations with time-dependent coefficients. Three criteria are introduced for model selection and threshold setting to identify the correct terms: the root mean square, total error bar, and group error bar. Moreover, three noise filters are integrated with the robust Bayesian sparse learning algorithm for better results with larger noise. Numerical results demonstrate that our method is more robust than sequential grouped threshold ridge regression and group Lasso in noisy situations through three examples.