Python 78 66 Bytes

lambda n:[i*[0]+[n]+(n+~i)*[0]for i in range(n)]+[n*[1+(n+1)**.5]]

Surely can be improved, especially at handling n=1```. (How is that even a simplex?) Just realized that's not necessary. Can probably be improved still ^^

Try it online!

[i*[0]+[1]+(n+~i)*[0]for i in range(n)] creates identity matrix. All points have distance sqrt(2) from each other. (thanks to Rod for improving)

Now we need a n+1-th point with the same distance to all other points. We have to choose (x, x, ... x).

Distance from (1, 0, ... ) to (x, x, ... x) is sqrt((x-1)²+x²+...+x²). If we want an n dimensional simplex this turns out to be sqrt((x-1)²+(n-1)x²), as we have one 1 and n-1 0s in the first point. Simplify a bit: sqrt(x²-2x+1+(n-1)x²) = sqrt(nx²-2x+1)

We want this distance to be sqrt(2).

sqrt(2) = sqrt(nx²-2x+1)
2 = nx²-2x+1
0 = nx²-2x-1
0 = x²-2/n*x+1/n

Solving this quadratic equation (one solution, other one works fine, too):

x = 1/n+sqrt(1/n²+1/n) = 1/n+sqrt((n+1)/n²) = 1/n+sqrt(n+1)/n = (1+sqrt(n+1))/n

Put that in a list n times, put that list in a list and join with identity matrix.

-4 Bytes thanks to Alex Varga:

Multiply each vector by n. This changes the creation of the identity matrix to lambda n:[i*[0]+[n]+(n+~i)*[0] (same length) and gets rid of the division by n in the additional point, so it becomes n*[1+(n+1)**.5], saving two brackets and the /n.

APL (Dyalog), 20 18 bytes

1 byte thanks to @ngn

∘.=⍨∘⍳⍪1÷¯1+4○*∘.5

Try it online!

R, 38 bytes

function(n)rbind(diag(n,n),1+(n+1)^.5)

Try it online!

Pari/GP, 37 bytes

n->concat((n+1)^.5-matid(n),1^[1..n])

Try it online!

Octave, 31 bytes

Saved 2 bytes thanks to Luis Mendo.

@(n)[n*eye(n);~~(1:n)+(n+1)^.5]

Try it online!

Ruby, 55 bytes

rather than returning similar magnitudes for all dimensions and using the formula (1+(n+1)**0.5)/n I scale up by a factor of n to simplify the formula to (1+(n+1)**0.5)

->n{(a=[n]+[0]*~-n).map{a=a.rotate}+[[1+(n+1)**0.5]*n]}

Try it online!

ungolfed in test program

A lambda function taking n as an argument and returning an array of arrays.

f=->n{
  (a=[n]+[0]*~-n).map{        #setup an array `a` containing `n` and `n-1` zeros. perform `n` iterations (which happens to be the the size of the array.)
  a=a.rotate}+                #in each iteration rotate `a` and return an array of all possible rotations of array `a`     
  [[1+(n+1)**0.5]*n]          #concatenate an array of n copies of 1+(n+1)**0.5
}

p f[3]                        # call for n=3 and print output

output

[[0, 0, 3], [0, 3, 0], [3, 0, 0], [3.0, 3.0, 3.0]]

Wolfram Language (Mathematica), 205 bytes

f1 = Sqrt[# (# + 1)/2]/# /(# + 1) & ;
f2 = Sqrt[# (# + 1)/2]/# & ;
simplex[k_] := {ConstantArray[0, k]}~Join~Table[
   Table[f1[n], {n, 1, n - 1}]~Join~{f2[n]}~Join~
    ConstantArray[0, k - n],
   {n, k}]

Simplex function in Mathematica Starting from {0,0,...]},{1,0,0,...]}, Placing first point at origin, Second point on x axis Third point in x,y plane, Fourth point in x,y,z space, etc. This progression reuses all the previous points, adding one new point at a time in new dimension

simplex[6]={{0, 0, 0, 0, 0, 0}, {1, 0, 0, 0, 0, 0}, {1/2, Sqrt[3]/2, 0, 0, 0, 
  0}, {1/2, 1/(2 Sqrt[3]), Sqrt[2/3], 0, 0, 0}, {1/2, 1/(2 Sqrt[3]), 
  1/(2 Sqrt[6]), Sqrt[5/2]/2, 0, 0}, {1/2, 1/(2 Sqrt[3]), 1/(
  2 Sqrt[6]), 1/(2 Sqrt[10]), Sqrt[3/5], 0}, {1/2, 1/(2 Sqrt[3]), 1/(
  2 Sqrt[6]), 1/(2 Sqrt[10]), 1/(2 Sqrt[15]), Sqrt[7/3]/2}}

Verification

In[64]:= EuclideanDistance[simplex[10][[#[[1]]]],simplex[10][[#[[2]]]]] & /@ Permutations[Range[10],{2}]//Simplify
Out[64]= {1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1}

Jelly, 11 bytes

‘½‘÷ẋW
=þ;Ç

Try it online!

Works by generating the identity matrix of size N and concatenating it with the list generated by repeating N times the singleton √(N + 1) + 1, divided by N.

‘½‘÷ẋW – Helper link (monadic). I'll call the argument N.

‘      – Increment N (N + 1).
 ½     – Square root.
  ‘    – Increment (√(N + 1) + 1).
   ÷   – Divide by N.
    ẋ  – Repeat this singleton list N times.
     W – And wrap that into another list.

––––––––––––––––––––––––––––––––––––––––––

=þ;Ç   – Main link.

=þ     – Outer product of equality.
  ;Ç   – Concatenate with the result given by the helper link applied to the input.

JavaScript (ES7), 70 bytes

n=>[a=Array(n++).fill((1+n**.5)/--n),...a.map((_,i)=>a.map(_=>+!i--))]

Port of @PattuX's Python answer.

Wolfram Language (Mathematica), 46 bytes

IdentityMatrix@#~Join~{Table[1-(#+1)^.5,#]/#}&

Try it online!