Maximal 2-distance Sets
Mathematica, 220 218 bytes

For[n=1,1&gt;0,n++,For[m=1,1&gt;0,m++,x=Array[a,{m,n}];d=Table[(Norm[x[[i]]-x[[j]]]^2-1)*(Norm[x[[i]]-x[[j]]]^2-y)==0,{i,1,m-1},{j,i+1,m}];If[Not[CylindricalDecomposition[Flatten@{d,y&gt;0},Flatten@{y,x}]],Break[]]];Print[m-1]]


Try it online! Un-golfed version:

For[ n = 1, 1&gt;0, n++,
  For[ m = 1, 1&gt;0, m++,
   x = Array[a, {m, n}];
   d = Table[(Norm[x[[i]] - x[[j]]]^2 - 1)*(Norm[x[[i]] - x[[j]]]^2 - 
         y) == 0, {i, 1, m - 1}, {j, i + 1, m}];
   If[Not[
      CylindricalDecomposition[Flatten@{d, y &gt; 0}, Flatten@{y, x}]],
      Break[]]
  ];
  Print[m - 1]
]


In this code, we apply a key tool from real algebraic geometry. A semialgebraic set is a subset of \$\mathbb{R}^d\$ defined by polynomial equalities and inequalities. We may identify the set of 2-distance subsets of \$\mathbb{R}^n\$ of cardinality \$m\$ with a semialgebraic subset of \$\mathbb{R}^{mn+2}\$. We force one of the distances to equal 1 without loss of generality, which means our semialgebraic set is a subset of \$\mathbb{R}^{mn+1}\$. We denote the other distance (squared) by \$y\$.

For which values \$Y\subseteq\mathbb{R}\$ of \$y\$ does there exists a size-\$m\$ subset of \$\mathbb{R}^n\$ with squared distances \$\{1,y\}\$? This can be answered with the help of Tarski--Seidenberg, which says that the projection of any semialgebraic set is semialgebraic. In our case, this means there is some univariate polynomial \$p\$ and relation \$*\in\{=,&lt;,\leq\}\$ such that \$Y=\{y\in\mathbb{R}:p(y)*0\}\$. This projection may be accomplished using cylindrical algebraic decomposition (CAD), which enjoys an implementation in Mathematica.

Given a set of polynomial equalities and inequalities, as well as an ordered list of variables, Mathematica's implementation returns either an explicit description of every point in the corresponding semialgebraic set, or False if the set is empty. In the above code, we run CAD for increasingly larger values of \$m\$ until we get False, at which point we know that \$s_n=m-1\$. Sadly, the runtime of CAD is doubly exponential in the number of variables, and so we shouldn't expect this approach to determine \$s_9\$ any time soon. In fact, a minute of runtime only gives \$s_1=3\$.