Tl;dr: The answer to your question is, Yes, but you picked the wrong weights on the graph, which is why the eigenvalues were off.
Before I get into some generalities about how to think about this kind of approximation problem, in addition to the paper Steve Huntsman posted above, you may find this paper by Burago-Ivanov-Kurylev interesting. The main theorem on eigenvalue approximation is Theorem 1:
For every integer $n\geq 1$ there exist positive constants $C_n$ and $c_n$ such that the following holds. Let $M$ be an $n$-dimensional Riemannian manifold with sectional curvature absolutely bounded by $K$, diameter bounded by $D$, and injectivity radius bounded below by $i_0$. Let $\Gamma(X,\mu,\rho)$ be a weighted graph defined by taking an $\epsilon$-net $X$, connecting all vertices that are less than $\rho$ apart, and weighting the edges and vertices to approximate the volume form on $M$. (The procedure is described starting at the bottom of page $3$.)
Suppose $\rho < i_0/2$, $K\rho^2 < c_n$, and $\epsilon/\rho < \min\{1/n,1/3\}.$ Then for every positive integer $k$, such that $\rho\lambda_k(M) < c_n$, one has
$$ \big|\lambda_k(\Gamma) - \lambda_k(M)\big| \leq C(n,K,D,i_0,k)(\epsilon/\rho + \rho) $$
where $C$ depends on the index $k$, the dimension $n$, and the bounds on curvature, diameter, and injectivity radius.
This theorem codifies your intuition of choosing points, connecting them at small scales, and then using the graph Laplacian to estimate the Laplacian on the manifold. However, we need to do better than use the adjacency matrix, since simply recording "these points are nearby" doesn't reflect the geometry of the manifold well enough to pick up on the eigenvalues of the Laplacian. So Burago, Ivanov, and Kurylev introduce a measure on the vertices of the graph and a weighting on each edge in order to make sure that the $L^2$ inner product and energy inner product on the graph correctly approximate the $L^2$ inner product on $M$ and the Laplacian on $M$.
More generally, this question about discretizing the Laplace eigenvalue problem is really a subset of numerically approximating solutions to PDEs. In that vein, you may also be interested in finite element approximation and the method of particular solutions. For finite elements, this is a reference, and there are good books by Brenner-Scott and Strang-Fix. For software, you could check out deal.ii (C++) and pydec (python, can be used to find eigenforms as well as eigenfunctions). For the method of particular solutions, perhaps consult this paper by Betcke-Trevethen, which has been adapted to hyperbolic surfaces in this paper by Strohmaier-Uski. Their software is also available on Alex Strohmaier's website (it's in Fortran).
Let me take a little bit to describe what all of these approaches have in common. First, let's set up the functional analysis. Recall that on a smooth compact Riemannian manifold $M$ the Laplace operator on smooth functions is defined as $\Delta = -\operatorname{div}\operatorname{grad}$.
If $M$ is closed, then $\Delta$ is essentially self-adjoint on smooth functions. If $M$ has boundary, then by restricting to certain smooth functions that satisfy self-adjoint boundary conditions (Dirichlet or Neumann, say), we produce a domain on which $\Delta$ is self-adjoint. In any of these cases, call the domain $\mathcal{D}$. Notice that $\mathcal{D}$ is a subset of $L^2(M)$. Denote the $L^2$ inner product by $n(\cdot,\cdot)$.
Integrating $\Delta$ by parts for a function in $\mathcal{D}$ produces an energy functional $a$, defined by $$a(v,w) = n(v,\Delta w) = \int_M \langle \nabla v,\nabla w\rangle\ dV.$$
Complete the domain $\mathcal{D}$ with respect to the norm $a(v,v) + n(v,v)$. In the case of the Laplacian, this is the $1$-Sobolev norm. Call the resulting space $H$. (When done in general for a self-adjoint operator, this process is called the Friedrichs extension.)
So far, so good. Eigenvectors and eigenfunctions of the Laplacian then have the following property: $v$ is a $\lambda$-eigenfunction if and only if, for all $w\in H$, we have $a(v,w) = \lambda n(v,w)$. (See this by integrating by parts.)
The intuitive goal of discretization for the eigenvalue problem is,
Find a finite-dimensional subspace $V$ of $H$ so that the eigenvalue problem restricted to $V$ closely approximates the eigenvalue problem on $H$.
By doing so, one reduces the infinite-dimensional eigenvalue problem to a finite-dimensional computation.
As an example, I'll describe the basic construction of the finite element approximation of the Neumann Laplacian on a bounded domain $\Omega$ in $\Bbb{R}^n$. (I believe this is called, more generally, the conforming Galerkin method.) Take some triangulation $\Gamma_h$ of $\Omega$ such that the diameter of every triangle in $\Gamma_h$ is less than $h$, and such that the smallest angle is no less than, say, $\pi/10$. Label the vertices of this triangulation $x_i$, for $i=1,\ldots,N$. For each $i$, define the piecewise-linear function $\phi_i$ to be the unique piecewise linear functions on $\Omega$ which has the property that $\phi_i(x_j) = \delta_{ij}$.
Let $S^h$ be the subspace of $H$ spanned by the $\phi_i$. The restricted eigenvalue problem is this: find $u^h$ and $\lambda^h$ such that $a(u^h,v^h) = \lambda^h n(u^h,v^h)$ for all $v^h\in S^h$. One can show (and indeed, Strang and Fix do this fairly explicitly in their book) that the relative error of the approximations that result is controlled by a power of $h$:
$$ \frac{|\lambda - \lambda^h|}{\lambda} \leq Ch^r $$
where $r$ is related to the dimension $n$ and the order of the operator (in this case, two), and $C$ depends on quantities like the ellipticity and coercivity of $\Delta$ and the minimum angle of the mesh (this is why I disallowed "bad triangles" by bounded minimum angles away from zero).
To actually perform computations, one solves for the matrices $A$ and $N$ by letting $A_{ij} = a(\phi_i,\phi_j)$ and $N_{ij} = n(\phi_i,\phi_j)$. Then one solves the finite-dimensional eigenvalue problem $Au = \lambda Nu$. (Note that the approximate eigenfunction is given by $u = \sum u_i\phi_i$.) This can be thought of as a weighted graph Laplacian and a weighted $L^2$ space on the graph $\Gamma^h$.
Burago-Ivanov-Kurylev do something similar, albeit they take more care in defining maps between $L^2$ of their graph and the domain of the Laplacian. In the method of particular solutions, one also defines an appropriate subspace, albeit in terms of a boundary discretization, rather than an interior discretization.
