The classical Gompertz model for tumor growth is given by the equations
$$\frac {dN}{ dt }=G (t) N(t)$$
$$\frac {dG}{ dt }=- \alpha G (t)$$
with initial conditions $N(0)=N_0$ and $G(0)=G_0$. Here $ N(t)$ is the number of cells at time $t $, $ G(t)$ is the "per capita" growth rate (that is to say, birth rate minus death rate at time $t$), and $\alpha $ is a positive constant. The model assumes that the tumor growth rate decays exponentially over time at rate $\alpha $, as a result of increasing death rates, decreasing reproduction rates, or both. Integrating the second equation leads to $G (t)=G_0 e^{-\alpha t} $, so that the model is usually written as
$$\frac {dN}{ dt }=G_0 \, e^{-\alpha t} N (t) $$
where the solution of this nonautonomous ODE is
$$ N (t)= N_0 \, e^{\alpha^{-1}(G_0-G(t))}$$
For practical purposes, however, the initial model is often simplified by setting $G(t) = (\beta− \alpha \log {N})$, where $\beta$ is an additional positive constant. This allows to write the following common formulation of the Gompertz model:
$$\frac {dN}{ dt }=N(\beta− \alpha \log {N})$$
An advantage of this simplification is that we have an autonomous equation, where the equilibrium points corresponds to the zeroes of the RHS. Thus, we have to solve
$$N(\beta− \alpha \log {N})=0$$
The solutions are $N=0$ and $N=e^{\beta/\alpha} $, which are then the equilibrium points (note that the equilibrium at $N=0$ is only theoretical, since in this point the equation is not defined). To assess the stability of these equilibrium points, considering that by definition $N \geq 0$, we can note that
$$N(\beta− \alpha \log {N}) >0 \,\,\,\text{for}\,\,\, 0<N <e^{\beta/\alpha}$$
$$N(\beta− \alpha \log {N}) <0 \,\,\,\text{for}\,\,\, N >e^{\beta/\alpha}$$
So, if $$ 0<N <e^{\beta/\alpha}$$ then $N' >0$, which means that the solutions increase towards the equilibrium point $N=e^{\beta/\alpha} $. Also, if $$N >e^{\beta/\alpha}$$ then $N'<0$, which means that the solutions decrease towards the equilibrium point $N=e^{\beta/\alpha} $. Because all solutions near $N=0$ increase and then move away from this point, and because all solutions near $N <e^{\beta/\alpha}$ converge towards this point, we can conclude that the theoretical equilibrium at $N=0$ is unstable, whereas the equilibrium at $N=e^{\beta/\alpha} $ is stable.
Regarding the simpler von Bertalanffy model, its most classic formulation is
$$\frac {dM}{dt}= \alpha M^{2/3}−\beta M$$
where $M $ is the tumoral mass, and $\alpha $ and $\beta $ can be interpreted as per-capita birth and death rates, respectively. The model assumes that the birth rate is proportional to the surface of the tumor, while the death rate is proportional to the mass of the tumor. Typically, this model is simplified assuming that the growth is homogeneous in the three dimensions, so that calling $L $ a linear dimension of the tumor we can write $$M (t)=k L (t)^3$$ where $k $ is a constant. From this, it follows $$dM=3k L^2 dL $$ and substituting in the original equation we get
$$\frac {dL}{dt}= \frac {1}{3}(\tilde {\alpha} -\beta L)$$
where $ \tilde {\alpha}= \alpha/k^{1/3} $. This shows that
$$L = \frac {\tilde {\alpha}}{ \beta} $$
is the only equilibrium point of the von Bertalanffy model, and that it is asymptotically stable.
Best Answer
You can derive their analytic expression by first taking the growth rate as a composite function of time, $(g\circ x)(t)$, and using this to derive the tumor cell number function, $x(t)$. I'm using this notation as opposed to $g(x)$ to emphasise that $x$ is also a function of time.
$$\begin{align}\frac{\mathrm{d}x(t)}{\mathrm{d}t}&=(g\circ x)(t) \\ \end{align}$$
Using the assumption that the growth rate is Gompertzian, $(g\circ x)(t)=ax(t)\ln\left(\frac{b}{x(t)}\right)$, with constants, $a,b\in\mathbb{R}$, we can plug this into the previous equation to get $\displaystyle \frac{\mathrm{d}x(t)}{\mathrm{d}t}=ax(t)\ln\left(\frac{b}{x(t)}\right)$. Shift all $x(t)$ terms to one side (LHS) and all $t$ terms to the other. Then integrate wrt. $t$ and use the fundamental theorem of calculus to deal with the derivative term.
$$ \begin{align}&\frac1{x(t)\ln\left(\frac{b}{x(t)}\right)}\frac{\mathrm{d}x(t)}{\mathrm{d}t}=a \\ &\int \frac1{x(t)\ln\left(\frac{b}{x(t)}\right)}\frac{\mathrm{d}x(t)}{\mathrm{d}t}\,\mathrm{d}t=\int a\,\mathrm{d}t \\ &\int \frac1{x(t)\ln\left(\frac{b}{x(t)}\right)}\,\mathrm{d}x(t)=at+C \end{align}$$
Then, since $\displaystyle \int \frac{1}{u\ln(b/u)}\,\mathrm{d}u=-\ln\ln\frac{b}{u}+C$, we have $\displaystyle -\ln\ln\frac{b}{x(t)}=at+C$ and therefore,
$$ x(t)=b\exp\left(-\exp\left(-\left(at+C\right)\right)\right)$$
setting $x(0)=N_0$ and solving for $C=-\ln\ln\left(\frac{b}{N_{0}}\right)$ shows that this agrees with the solution of the paper and intuitively makes sense since we'd expect the curve to be sigmoidal. The paper was slightly disingenuous when it suggested that you integrate wrt. $x$, since there was some preceding algebra that you'd have needed to do, and it's sort of an integration wrt. $t$ in disguise. Your solution of $a\left(\frac{\ln\left(\frac{b}{x}\right)x^{2}}{2}+\frac{x^{2}}{4}\right)$ suggests that you were integrating $g$ outright wrt. $x$ without taking this into account. If you do this, you get an expression for $\displaystyle \int \frac{\mathrm{d}x(t)}{\mathrm{d}t}\,\mathrm{d}x $, which doesn't help much since we essentially want the differentials to cancel. We can also see that it was intuitively wrong, since it exploded to $-\infty$ with large $t$.
To approach the problem numerically with Euler's method, we have a system defined by $\dfrac{\mathrm{d}x}{\mathrm{d}t}=g(x),\ x(0)=N_0$ with the Gompertz function mentioned earlier. This can be given an approximate solution $\big((0,N_0),(h,y_1),(2h,y_2),\ldots,(nh,y_n)\big)$ where $y_{n+1}=y_n+hg(y_n)$ and $y_0=N_0$. The error of the function at each point could be defined as $\varepsilon=|y_i-x(ih)|$ using the known solution $x(t)$. The figure below shows the approximate solution you'd get for the curve defined by the parameters $a=1,b=2,N_0=0.1$, using points of $h=1.5$.