In our calculation of Sx(f) from Umax prior to this, we were missing a factor of 4. The correct formula from Liu and Thorne's formula 58 is Sx[f] := 4*Kb*T*Umax*Phi/Pi*f*F0^2). When we correct this formula, we get results which agree fairly well with the analytic results. The Mathematica script attached was used to perform both calculations of the thermal noise spectrum (using SxComsol and SxFTMThorne).
SxCOMSOL(100 Hz) = 7.99735*10^40 m^2/Hz
SxFTMThorne(100 Hz) = 7.80081*10^40 m^2/Hz
These results differ by about 2.5%. We also need to verify that the result converges for increasing mesh size.
Number of Elements

Umax (*10^10 J)

1224

1.48859

7011 (used above)

1.51589

13306

1.51891

I attempted to run the simulation with 24360 elements in the mesh, but the computer I was running it on repeatedly crashed.
edit: Running on a more powerful computer I got the following additional values.
Number of Elements

Umax (*10^10 J)

1224

1.48859

7011 (used above)

1.51589

13306

1.51891

24181

1.52169

61772

1.52281

Given these additional values, it appears that our simulation will converge to a value for increasing mesh size and that it agrees fairly well with the analytic result. 