Probabilistic numerics methods are a novel approach to quantifying the approximation errors in numerical computations as probabilistic uncertainties. A recent method that was developed is the Bayesian Finite Element Method, which aims to determine the discretization errors along a coarse mesh probabilistically. This work analyzes the use of priors in this method. It is shown that the priors in the right-hand side or forcing term of the partial differential equation are superior to applying them directly in the solution. It is demonstrated that the maximum log-likelihood is not an appropriate estimator of hyperparameters, and as solutions, hyperparameters are optimized using objective functions to capture the discretization error with the posterior deviation. The advantages of non-stationary priors are studied in order to have standard deviations that show the error along the mesh. Additionally, this work examines how the optimal hyperparameters change for mesh refinements and different arrangements of elements in the same problem. Finally, the research delves into determining the approximate number of samples necessary for ensembling covariance matrices and obtaining similar results.