Based on the comments below your post you may reach the conclusion that # DOFs and speed have no correlation whatsoever - this is not true. Keeping all other things fixed and increasing the number of unknowns it is clear that you would need to access more memory and perform more arithmetic operations. You don't need to trust me on that - you can test linear vs quadratic triangular elements where you have fixed the number of elements to be the same. The latter results in a slightly denser matrix and has more unknowns (# of vertices + # of edges versus just # of vertices for linear).
However the reason for choosing higher order polynomials is to get a better decrease rate of the error with the decrease of the size of the elements, so the goal for choosing higher order elements is slightly different. A fairer comparison there would be to compare linear elements on a mesh with $n$ vertices vs quadratic element on a mesh where the number of vertices plus edges sums up to $n$. Then you would expect the latter to get closer to the true solution as you decrease the element size.
In FEM if you have $m$ Dirichlet nodes and $(n-m)$ independent nodes, you get $n-m$ unknowns, and a (sparse) matrix that is $(n-m)\times (n-m)$. While you can grow the number of unknowns in order to get a better approximation, this generally also makes the solving process slower. If $m\ll n$ then you may opt to use Green's functions for the solution instead, which result in systems of size $m\times m$. These systems are typically dense, but as long as $m$ is small enough you can use a direct solver on those. The main practical issue there is whether you can construct the Green's functions efficiently. For example one can do so for the Poisson equation on a rectangular domain with a regular grid and with zero Dirichlet, reflecting, or periodic boundary conditions at the rectangle's boundaries, by using respectively the discrete sine, cosine, or Fourier transform (then the resulting solution strategy is typically termed fast Poisson solver. In the general case one has to compute columns of the generalised inverse of the non-reduced FDM/FEM/FVM matrix in order to get the Green's functions, so this can be more expensive than just running an iterative solver on the sparse formulation.
Discretisations based on the Green's function formulation are typically termed BEM or RBF based solutions in case you want to look those up. Typically a closed form solution for the Green's functions is not known for arbitrary meshings over arbitrary geometry however.
