A decent part of this answer is based on the introduction to Kroupa & Weidner (2005), though I’ve obviously gone into a lot more depth on all of the references.
Our story starts, as do many concerning stellar astrophysics, with Sir Arthur Eddington. In his 1926 book, The Internal Constitution of the Stars, he derived the Eddington luminosity, the maximum luminosity $L$ a star of mass $M$ can reach (Chapter 6, pages 114-115). His derivation goes along the following lines:
I. Take the equation of hydrostatic equilibrium and the equation of radiative equilibrium:
$$\frac{\mathrm{d}P}{\mathrm{d}r}=-g\rho\tag{1a}$$
$$\frac{\mathrm{d}p_R}{\mathrm{d}r}=-\frac{k\rho H}{c}\tag{1b}$$
The relevant variables are pressure ($P$), radius ($r$), gravitational acceleration ($g$), density ($\rho$), radiation pressure ($p_R$), mass coefficient of absorption ($k$), radiative flux per time ($H$), and the speed of light ($c$). Combining $(1\text{a})$ and $(1\text{b})$ yields
$$\mathrm{d}p_R=\frac{kH}{cg}\mathrm{d}P\tag{1c}$$
II. At some radius $r$, the luminosity $L_r$ and enclosed mass $M_r$ can be related by
$$\frac{L_r}{M_r}=\frac{\eta L}{M}\tag{2a}$$
where $L$ and $M$ are the luminosity and enclosed mass at the radius of the star, and $\eta$ is some function of $r$, increasing inward from $\eta(R)=1$ at the stellar radius $R$. Given that
$$H=\frac{L_r}{4\pi r^2}\tag{2b}$$
$$g=\frac{GM_r}{r^2}\tag{2c}$$
we have
$$\frac{H}{g}=\frac{L_r}{4\pi GM_r}\tag{2d}$$
Putting this back into $(1\text{c})$, we find
$$\mathrm{d}p_R=\frac{L\eta k}{4\pi cGM}\mathrm{d}P\tag{2e}$$
III. As temperature and density increase towards the center of the star, so does the pressure due to the matter, $p_G$. Therefore, $\mathrm{d}p_G>0$. Furthermore, given that $P=p_G+p_R$, $\mathrm{d}p_R<\mathrm{d}P$. This means that $(2\text{e})$ yields
$$\frac{L\eta k}{4\pi cGM}<1\tag{3}$$
which is the criterion leading to the Eddington luminosity. There are, of course, other ways of getting this criterion, but I thought I’d give Eddington’s original one, in all its mathematical glory.
Using a suitable mass-luminosity relation for massive stars, we can then establish the mass of a star at the Eddington limit. Eddington himself took it to be in the range of 60-70 solar masses ($M_{\odot}$), though today a value of somewhere around 120 solar masses is more appropriate.
Let us take a detour to a lesser-known figure, Paul Ledoux. In 1941, Ledoux analyzed vibration modes in stars due to the usual perturbations in density, pressure, radius, temperature, etc. He came up with the stability condition of
\begin{align}
A_k&=\!
\begin{aligned}[t]
\int_0^M\frac{\delta \rho_k}{\rho}\biggl[(\Gamma_3-1)\delta_k\left\{\epsilon_1+\epsilon_2-\epsilon_3-\frac{\mathrm{d}}{\mathrm{d}m}[4\pi r^2(F_1+F_3)]\right\}&\\
-\frac{2}{3}\delta_k\left[4\pi r^2\bar{C}\frac{\mathrm{d}P}{\mathrm{d}m}+\epsilon_2+\frac{\mathrm{d}}{\mathrm{d}m}[4\pi r^2F_2]\right]\biggr] \mathrm{d}m&<1\\
\end{aligned}\\
\end{align}
for the $k$th mode of vibration. I’m not going to explain all the variables because that’s not quite important; the important takeaway is that Ledoux took turbulent pulsations into account. His conclusion is that an exact model “probably” would lead to a limit of about 100 solar masses; using certain inexact assumptions, he found a limit of 128 solar masses.
Ledoux’s analysis laid the groundwork for the work of Schwarzschild & Härm (1958). Their stability criterion isn’t necessarily simpler, but it can be written more compactly. Specifically, the coefficient of stability, $K$, defined as
$$K=\frac{1}{2}\frac{L_P}{E_P}$$
must be negative to ensure stability against pulsations. A positive $K$ means that the pulsation amplitude is increasing; a negative $K$ means that the pulsation amplitude is decreasing.
$E_P$ is the pulsation’s energy, while $L_P$ is the rate of gain of pulsation energy and can be expanded as
$$L_P=\overbrace{L_{PN}}^{\text{nuclear}}-\overbrace{L_{PH}}^{\text{heat leakage}}-\overbrace{L_{PS}}^{\text{progressive waves}}$$
Here, $L_{PN}$ represents the rate of energy gained, while $L_{PH}$ and $L_{PS}$ represent the rate of energy lost. All of the above quantities can be calculated through some relatively simple expressions (see Equations 9-12 and 15-22). The upshot of all this is that $K$ becomes negative at birth for stars greater than 60 solar masses. This can be figured out by writing $L_P$ and $E_P$ as functions of mass, $M$, and age, $\tau$.
Now, interestingly enough, the critical age ($\tau_{cr}$) can be written as a function of mass:
$$\tau_{cr}=0.05\left(\frac{M}{M_{\odot}}-60\right)$$
where $\tau_{cr}$ is in millions of years. This means that a star of, say, 62 solar masses (to take the authors’ example) will evolve to a stable state in a quarter of a million years. We can also determine whether or not, in this time, the star’s instability will become too great and destroy it. It turns out that this is the case for stars with masses greater than 65 solar masses - putting the upper limit for the mass of a star at 65 solar masses.
Here is a graphical representation from their paper, Figure 1:

Even later work on the same topic was done by Ziebarth (1970), among others, who extended the models to study different metallicities and compositions (Schwarzschild & Härm) focused largely on stars with compositions similar to that of the Sun). His computations found a wide range of upper mass limits - 10 solar masses for pure helium stars, and 200 solar masses for pure hydrogen stars. Most stars fall in the middle, and so will have different limits.
The actual formation of massive stars also puts constraints on the mass. Kroupa & Weidner mention Kahn (1974), who studied how radiation pressure from a protostar could drastically lower accretion rates, stopping the star from continuing to significantly grow. As applied to a young Population I star, his simplest model comes out to a limit of about 80 solar masses, although different models of the “cocoon” yield different results.
I’ll add one final note on theory. Population III stars, the hypothetical first stars in the universe, are expected to have been extremely massive; as such, they would be excellent candidates for testing the upper mass limits. According to simulations by Hosokawa et al. (2011), mechanisms similar to those discussed by Kahn would have stopped accretion at stellar masses around 43 solar masses - a surprisingly low figure, given the expectations of how massive Population III stars should be. Additionally, as argued by Turk et al. (2009), sufficiently massive stars could fragment; in the case studied, a 50 solar mass star broke apart into two smaller core fragments.
Something I realized just now, a couple months after writing this, is that all of this assumes that the star is spherically symmetric. Most stellar models involve spherically symmetric, non-rotating stars, which allows us to make some assumptions such that the equations of stellar structure depend purely on $r$, the radial coordinate.
We have, however, seen stars - not stellar just simply stellar remnants like pulsars, but even main-sequence stars - that rotate rapidly and are thus non-spherical. Vega, for instance, has an equatorial radius 19% larger than its polar radius. If a star of mass $M$ is rotating, the equations of stellar structure should be different, and so some of the above results should also be different. I'm not sure how important this is for various theoretical limits.