Consider a wire of length $1$cm with constant density $\mu$(g/cm), containing two point-like masses $m_1$ and $m_2$ grams located at $s_1=\frac{1}{3}$ and $s_2=\frac{2}{3}$, assuming the wire is modeled as a real line [0,1]. Find a function $\rho:[0,1]\to \mathbb{R}$ such that a total mass of any portion $(a,b]\subset [0,1]$ of the wire is given by $M(a,b)=\int_a^b d\rho$.
By this setting, our goal is to find a function $\rho$ so that
$M(a,b)=\begin{cases} \mu(b-a) & s_1,s_2\notin (a,b] \\ \mu(b-a)+m_1 & s_1\in (a,b], s_2\notin (a,b] \\ \mu(b-a)+m_2 & s_1\notin (a,b], s_2\in (a,b] \\ \mu(b-a)+m_1+m_2 & s_1,s_2\in (a,b] \end{cases}$
If there are no point masses, then I think $\rho(x)=\mu x$ is our desired function. Also, this $\rho$ should be a monotonically increasing function since this integral is Riemann-Stieltjes. Because I've never encountered this type of modeling problem, I got stuck here.
EDIT:
I naively guess
$\rho(x)=\begin{cases} \mu x & 0\leq x<\frac{1}{3} \\ \mu x+m_1 &\frac{1}{3}\leq x<\frac{2}{3} \\ \mu x +m_1+m_2 & \frac{2}{3}\leq x \leq 1 \end{cases}$
For a partition $P=${$x_0,...,x_n$} such that $x_{i-1}<s_1\leq x_i$, then $\Delta \rho_i=\rho(x_i)-\rho(x_{i-1})=\mu x_i+m_1-\mu x_{i-1}=m_1+\mu (x_i-x_{i-1})$. Letting $\Delta x_i \to 0$, then I can get the desired mass. What do you think about this?
To complement your approach. If you think about the fact that (for a continuously differentiable function $f$, one has $$ f(b) - f(a) = \int f'(x) dx, $$ then one should think about $d\rho$ as being something like $M'$. Don't take that literally, since $M$ is a function of two variables and not differentiable.
The next step is to note that only $M(0,b)$ matters, because $M(a,b) = M(0,b) - M(0,a)$ and the integral is linear. Now, $$ M(0,x) = \mu x + m_1\mathbb{1}_{x > x_1} + m_2\mathbb{1}_{x > x_2}. $$ This is already the function that you found.
We can put it in the form I wrote before by using $M(0,0) = 0$ so that we really want to compute $$ \rho(x) - \rho(0) = \int_0^x d\rho $$ for $\rho(x) := M(0,x)$. You might have seen that in a generalized sense, this makes sense with $$ \rho' = \mu\mathbb{1}_{[0,1]} + m_1\delta_{x_1} + m_2\delta_{x_2} $$ which is exactly the distribution that you want to see on the string.
Edit: Note that you don't need to start your integral from 0. You could have started with any another point, since that shifts $\rho$ only by a constant which we do not see in the differential.