<p>Approximate solutions of partial differential equations (PDEs) obtained by neural networks are highly affected by hyperparameter settings. For instance, the model training strongly depends on loss function design, including the choice of weight factors for different terms in the loss function, and the sampling set related to numerical integration; other hyperparameters, like the network architecture and the optimizer settings, also impact the model performance. On the other hand, suitable hyperparameter settings are known to be different for different model problems and currently no universal rule for the choice of hyperparameters is known. In this paper, for second order elliptic model problems in two to four spatial dimensions, various hyperparameter settings are tested numerically to provide a practical guide for efficient and accurate neural network approximation. While a full study of all possible hyperparameter settings is not possible, the focus is on the formulation of the PDE loss as well as the incorporation of the boundary conditions, the choice of collocation points associated with numerical integration schemes, and various approaches for dealing with loss imbalances on various model problems; in addition to various Poisson model problems, also a nonlinear and an eigenvalue problem are considered. The numerical results indicate that augmented Lagrangian balancing and more accurate Gaussian quadrature can substantially improve the performance of the deep Ritz method. Moreover, when applicable and suitably tuned, the deep Ritz method is often more accurate, stable, and computationally efficient than PINNs on more complex test cases, including three-dimensional, nonlinear, and eigenvalue problems. For the four-dimensional Poisson problems, however, while PINN-based approaches incur higher computational cost, there is no clear overall accuracy advantage for either approach.</p>