My primary experience is with crystal structures, and there is only a finite number of point symmetries that show up in a crystal. So, the algorithm I would use is slightly different than what you'd use in a molecule. But, it is unlikely with a large molecule that the continuous symmetries will show up, like the axial symmetry in H$_2$ or CO$_2$, so the methods should overlap quite well. When determining the symmetry in a system, there are two different, but related, symmetries to consider: local and global.
Local Symmetry
Local symmetry is the symmetry of the local environment around a specific point. In particular, the symmetry at each atomic location determines the local atomic splitting and to some extent the chemical environment, and is a subgroup of the global symmetry. For example, in benzene the local symmetry consists of two reflection planes and a $C_2$ axis ($180^\circ$ rotation symmetry). (Obviously, only two of the operations are necessary to generate the entire local point group.)
From an algorithmic perspective, what we've done is to first find the nearest neighbors of the target atom, and then enumerate all the ways we can rotate that environment about the central atom and have it remain the same. More mathematically, it is solving for all orthogonal matrices, $\mathbf{A}$, such that
$$\mathbf{A}(\vec{x}_i - \vec{x}_c) = \vec{x}_j - \vec{x}_c$$
where $\vec{x}_i$ and $\vec{x}_j$ are the positions of atoms of the same species and $\vec{x}_c$ is the position of the central, or target, atom. But, I'd look at simpler forms first, like whether or not a reflection plane exists, prior to trying to solve for $\mathbf{A}$ in general.
Another thought is to use the angular momentum matrices as generators of rotation, then
$$\mathbf{A} = \exp(i \phi \hat{n}\cdot\vec{\mathbf{L}})$$
where $\hat{n} \in \mathbb{R}^3$ is a unit vector about which a rotation of with angle $\phi$ is performed, and $\vec{\mathbf{L}} = (\mathbf{L}_x,\ \mathbf{L}_y,\ \mathbf{L}_z)$ is the vector of three dimensional angular momentum matrices. $\mathbf{A}$ would then have only 3 unknowns.
Global Symmetry
Where the local symmetry determines the environment around a single atom, the global symmetry dictates how the atoms interchange with each other. The first step in determining the global symmetry is to determine the equivalent atoms. First, determine the types of and the relative directions to the nearest neighbor (and second nearest, or higher, if desired) atoms. Two atoms are then equivalent, if their neighbors have the same spatial arrangements. This is straightforward to calculate.
The second step is roughly the same as that found in the local symmetry case, except that the center of mass of the molecule is likely the symmetry center. At this point, if the local symmetries have been determined, only a few unique operations may need to be found to generate the entire group. For example, in the B20 crystal structure, each atom has a local $C_3$ symmetry, and the full point group is generated by including a 2-fold ($180^\circ$ rotation) screw axis which transforms one atom into another. In benzene, two operations are required: a 6-fold ($60^\circ$) rotation through the central axis and a reflection plane bisecting a bond.
Edit: For the B20 structure, you can use two of the $C_3$ axes, instead, to generate the full group. This should allow you to avoid having to figure out a way to automatically determine the screw axis.
Caution: A caution on using the ideas in the local symmetry section in the global section, to be a symmetry operation, the environment must also be transformed. So, that if you find $\mathbf{A}$, from above, it will only give a candidate symmetry as the transformation may not similarly change the environment appropriately, and further checks are needed. For instance, if the benzene ring had hydrogen atoms sticking out of the plane of the ring along one side, then the reflection plane bisecting the carbon-carbon bond would be fine, but a $180^\circ$ rotation similarly bisecting the bond would not because it would not reproduce the local environment.
Edit - Translations: There is one other complication that the above discussion on local symmetry ignores: translations. Formally, the correct symmetry operation is
$$\mathbf{A}(\vec{x}_i - \vec{x}_c) + \vec{t} = \vec{x}_j - \vec{x}_c$$
where $\mathbf{A}$ and $\vec{x}_k$, as above, and $\vec{t}$ is an arbitrary translation. In a symmorphic crystal,
$$\vec{t} = n_1 \vec{a}_1 + n_2 \vec{a}_2 + n_3 \vec{a}_3$$
where $\vec{a}_i$ is a primitive lattice translation and $n_i \in \mathbb{Z}$, so the point group and translations are completely separable. In a non-symmorphic crystal, $\vec{t}$ can consist of non-primitive translations. The difference between the two is simply that for a symmorphic crystal, a single center of rotation can be found, but for non-symmorphic crystals, this is not true. A molecular system is likely to be "non-symmorphic" in this latter sense, and require the addition of translations to fully realize the group.
