Let $\mu$ be a Radon measure on $\mathbb{R}$. (i.e., a locally finite Borel measure)
Let $I$ be an interval of $\mathbb{R}$.
Let $\gamma : I \mapsto \mathbb{R}$ be a monotone function that is Lebesgue-Stieltjes integrable with respect to $\mu$: $$\int_I \gamma(x) \mathrm{d}\mu (x) < \infty.$$
I would like a simple argument to prove that there exists a sequence $(\gamma_n)_{n\in\mathbb{N}}$ of continuous and monotone functions $\gamma_n : I \mapsto \mathbb{R}$ such that $$\int_I \gamma_n(x) \mathrm{d}\mu (x) < \infty$$ and $$\lim_{n\to\infty} \int_I | \gamma(x) - \gamma_n (x)| \mathrm{d}\mu (x) =0.$$
I have thought I could use the Lusin theorem, but it would require to adapt its statement and proof to show the approximating functions are monotone. Besides, my hypotheses are much simpler ($\mathbb{R}$, no compacity of the support of approximating functions). So maybe a weaker theorem could give the answer? I am interested by any reference to such theorem.
N.B.: a proof of this is very common if $\mu$ is the Lebesgue measure, but here there is no absolute continuity assumption on $\mu$.
I don't think your version of Lusin's theorem can be used to give a short proof of the result, since we need some control over the continuous functions we get. Here's a solution, using a slightly weaker version of Lusin's theorem, and constructing the continuous functions more explicitly.
Let's assume $\gamma$ is non-decreasing (otherwise, work with $-\gamma$). First let's consider $I$ a compact interval of the form $I=[a,b]$, and let's deal with the general case later.
Let $\epsilon>0$. By Lusin's theorem, there exists a subset $A\subseteq I$ of measure $\geq \mu(I)-\epsilon$ such that $\gamma|_A$ (the restriction of $\gamma$ to $A$) is continuous. (This is the slightly weaker form of Lusin's theorem than the one you linked to). Moreover, by regularity of $\mu$, we can assume $A$ is closed.
We can also assume that $b\in A$: If this is not the case, then since $A$ is closed, $A\cup\{b\}$ is also closed and $\gamma|_{A\cup\{b\}}$ is continuous. Similarly, we can assume $a\in A$.
Now we will extend $\gamma|_A$ to a continuous, non-decreasing function on $I$. The complement $O=I\setminus A$ of $A$ is an open set, so it can be written as a countable union of disjoint open intervals $O=\bigcup_n (a_n,b_n)$. Note that $a_n,b_n\in A$ for all $n$.
Define $g=\gamma$ on $A$, and on each set $(a_n,b_n)$, define $g(x)=\gamma(a_n)+\frac{x-a_n}{b_n-a_n}(\gamma(b_n)-\gamma(a_n))$. Let's show that $g$ has the desired properties:
Suppose $x<y$ in $I$. We have some cases:
1.1. $x,y\in A$.
Then $g(x)=\gamma(x)<\gamma(y)=g(y)$.
1.2. $x\in A$, $y\not\in A$.
In this case, $y$ is in some set of the form $(a_n,b_n)$, and we must have $x\leq a_n$, so $$g(x)=\gamma(x)\leq\gamma(a_n)\leq g(y)$$
1.3. $x\not\in A$, $y\in A$.
This is similar to 1.2.
1.4. $x,y\not\in A$.
Then $x\in (a_n,b_n)$ and $y\in(a_m,b_m)$ for certain $n,m$. If $n=m$ then it is clear that $g(x)\leq g(y)$. If $n\neq m$, then the interval $(a_n,b_n)$ stands to the left of the interval $(a_m,b_m)$ (since they are disjoint and contain $x$ and $y$, respectively, and $x<y$), so $b_n\leq a_m$ and thus $$g(x)\leq g(b_n)=\gamma(b_n)\leq\gamma(a_m)=g(a_m)\leq g(y).$$
We leave the following result as an exercise:
Let's use this result: Let $x\in[a,b)$, and let's find a sequence as in the proposition. If $x\not\in A$, then $x\in(a_n,b_n)$ for some $n$ and it is easy enough to find such a sequence (e.g. $x_i=x+(b_n-x)/i$). So suppose $x\in A$. We have two cases:
2.1. There is a sequence $x_n\in A$ with $x_n>x$, $x_n\to x$.
Since $\gamma|_A=g|_A$ is continuous, we immediately have $g(x_n)\to g(x)$.
2.2. There is no sequence $x_n\in A$ with $x_n>x$, $x_n\to x$.
This means that there is some open interval of the form $(x,c)$ which does not intersect $A$, i.e., $(x,c)\subseteq I\setminus A$ and thus $(x,c)\subseteq(a_n,b_n)$ for some $n$. But $x\in A$, so $x=a_n$. Now it is again easy to construct the desired sequence (say $x_i=x+(b_i-a_i)/i$).
Therefore $g$ is right-continuous. Left continuity is similar and therefore $g$ is continuous.
Both $\gamma$ and $g$ are bounded on $I$ by $\max(|\gamma(a)|,\gamma(b)|)$, and they coincide on a set of measure $\geq\mu(I)-\epsilon$, thus $$\int_I|\gamma-g|d\mu\leq\max(|\gamma(a)|,\gamma(b)|)\epsilon$$ which is as small as necessary, by adjusting $\epsilon$.
Note that from this it follows that $\int_I|g|d\mu<\infty$.
The above solves the case of a compact interval $I$. If $I$ is of the form $I=[a,b)$, with $b\in(a,\infty]$, then we can write $I=\bigcup_{n=1}^\infty[a_n,a_{n+1}]$, where $a=a_1<a_2<\cdots$ and $a_n\to b$.
Given $\epsilon>0$, for each $n$ find a function $g_n$ on $[a_n,a_{n+1}]$, as above, which is continuous, monotone, $g_n(a_n)=\gamma(a_n)$, $g_n(a_{n+1})=\gamma(a_{n+1})$, and $\int_{[a_n,b_n]}|g_n-\gamma|d\mu<2^{-n}\epsilon$. Define $g$ on $I$ as $g|_{[a_n,a_{n+1}]}=g_n$, and we have the desired properties.
The last case is of intervals of the form $I=(a,b)$, with $-\infty\leq a<b\leq\infty$, and it follows from the previous case with arguments similar to the ones in the paragraph above.
Alternative solution:
We can assume that $I=[a,b]$ is compact and $\gamma$ is non-decreasing and non-constant. Moreover, up to normalization, we can assume $\gamma(a)=0$ and $\gamma(b)=1$. Let's start with the usual simple functions approximating $\gamma$: For each $n$ and $1\leq i<2^n$, take $A_{n,i}=\gamma^{-1}[\frac{i-1}{2^n},\frac{i}{2^n})$, and $A_{n,2^n}=\gamma^{-1}[\frac{2^n-1}{2^n},1]$.
Since $\gamma$ is monotone, each $A_{n,i}$ is an interval (possibly degenerate). Whenever $A_{n,i}$ has nonempty interior (i.e., it is not a point), take $a_{n,i}<b_{n,i}$ with $\operatorname{int}(A_{n,i})=(a_{n,i},b_{n,i})$. Notice that $a_{n,i}=0$ for precisely one $i$, and $b_{n,i}=1$ for precisely one $i$, and $\left\{a_{n,i}\right\}\setminus\left\{0\right\}=\left\{b_{n,i}\right\}\setminus\left\{1\right\}$.
Now take a decreasing sequence $\delta_n$ converging to $0$, in such a way that $a_{n,i}+\delta_n<b_{n,i}-\delta_n$ for all $n,i$. Let $B_{n,i}=(a_{n,i}+\delta_n,b_{n,i}-\delta_n)$.
Given $n$, define $g_n:I\to[0,1]$ as follows: $g_n(a_{n,i})=\gamma(a_{n,i})$, and $g_n(b_{n,i})=\gamma(b_{n,i})$ for all $i$. If $x\in B_{n,i}$ for some $i$, set $g_n(x)=i2^{-n}$. Then extend $g_n$ arbitrarily to some continuous, monotone function on $I$.
Let's show that $g_n$ converges pointwise to $\gamma$. For convenience, for each $n$ consider the partition $\mathscr{P}_n=\left\{A_{n,i}:i\right\}$. Note that $\mathscr{P}_m$ refines $\mathscr{P}_n$ whenever $m\geq n$.
Let $x\in I$ and $k\in\mathbb{N}$. We have two cases:
If $x=a_{k,i}$ for some $i$, then for every $m\geq k$ there will be some $j$ with $x=a_{m,j}$. This happens because $\mathscr{P}_m$ refines $\mathscr{P}_k$. So in this case we have $g_m(x)=\gamma(x)$ for all $m\geq k$. The same happens if $x=b_{k,i}$ for some $i$.
Now suppose $x$ is not any of the $a_{k,i}$ nor any of the $b_{k,i}$. Choose $i$ for which $x\in A_{k,i}$. Then $x$ lies in the interior of $A_{k,i}$, i.e., $x\in(a_{k,i},b_{k,i})$. Choose $N>k$ large enough such that $a_{k,i}+\delta_N<x<b_{k,i}-\delta_N$.
Now let's show that for any $m\geq N$, we have $|g_m(x)-\gamma(x)|\leq 2^{1-k}$. Again, we have some cases:
If $x$ is of the form $a_{m,j}$ or $b_{m,j}$ then $g_m(x)=\gamma(x)$, and we are done. So suppose $x$ is not of the form $a_{m,j}$ nor $b_{m,j}$. Again, there will be some $j$ for which $x\in\operatorname{int}(A_{m,j})=(a_{m,j},b_{m,j})$.
Take $p$ for which $b_{m,p}=b_{k,i}$. Again, since $\mathscr{P}_m$ refines $\mathscr{P}_k$, we have $A_{m,p}\subset A_{k,i}$, and the definition (and nonemptyness) of these sets garantees that $p2^{-m}\leq (i+1)2^{-k}$. Indeed, given $y\in A_{m,p}$, we have $$(p-1)2^{-m}\leq\gamma(y)\leq p2^{-m}\qquad\text{and}\qquad (i-1)2^{-k}\leq\gamma(y)\leq i2^{-k}$$ so $$p2^{-m}\leq\gamma(y)+2^{-m}\leq\gamma(y)+2^{-k}\leq (i+1)2^{-k}.$$
By our choice of $\delta$, we have $x<b_{n,i}-\delta_N\leq b_{m,p}-\delta_m$, so we can find some $y\in B_{m,p}$ with $x\leq y$. Since $g_m$ is monotone, $$g_m(x)\leq g_m(y)=p2^{-m}\leq (i+1)2^{-k}$$
With the same kind of reasoning (now working with the $a_{m,j}$), we can prove that $$(i-1)2^{-k}\leq g_m(x)$$ so we have $$(i-1)2^{-k}\leq g_m(x)\leq(i+1)2^{-k}$$ and also, now by definition of $A_{k,i}$, $$(i-1)2^{-k}\leq \gamma(x)\leq i2^{-k}\leq (i+1)2^{-k}$$ so $|g_m(x)-\gamma(x)|\leq 2^{1-k}$, as we wanted.
Therefore $g_m\to \gamma$ pointwise (everywhere!), and dominated convergence implies $\int g_m\to\int\gamma$. We did not assume anything about the measure (just that it was finite on $I$).