Sadly, I did not find any book in which specifically this result is
mentioned, so I will give a proof sketch below.
You might find that you already know some of the facts I use. In this case, skim. Otherwise, I have provided some links where you might find more information on the individual steps.
We start by approximating $f\in W^{1,1}\left(\left(0,1\right)\right)$
(the same works for any open interval of finite length) locally with
smooth functions.
Now, for $\varepsilon>0$, let $\varphi_{\varepsilon}:=\frac{1}{\varepsilon}\cdot\varphi\left(\frac{x}{\varepsilon}\right)$.
Then $\varphi_{\varepsilon}\in C_{c}^{\infty}\left(\left(-\varepsilon,\varepsilon\right)\right)$,
$\int_{\mathbb{R}}\varphi_{\varepsilon}\, dx=1$ and for every $g\in L^{1}\left(\mathbb{R}\right)$,
we have $\varphi_{\varepsilon}\ast g\xrightarrow[\varepsilon\downarrow0]{L^{1}\left(\mathbb{R}\right)}g$,
where the so-called convolution $F_{\varepsilon}:=\varphi_{\varepsilon}\ast g$
of $\varphi_{\varepsilon}$ and $g$ is given by
$$
\left(\varphi_{\varepsilon}\ast g\right)\left(x\right)=\int_{\mathbb{R}}\varphi_{\varepsilon}\left(x-y\right)g\left(y\right)\, dy.
$$
On convolution in general, see also here Why convolution regularize functions?
and here https://mathoverflow.net/questions/5892/what-is-convolution-intuitively
- To show this, we use that $\left\Vert g-L_{x}g\right\Vert _{1}\xrightarrow[x\rightarrow0]{}0$
for every $g\in L^{1}\left(\mathbb{R}\right)$, where $\left(L_{x}g\right)\left(y\right)=g\left(y-x\right)$.
This is shown for example here proof that translation of a function converges to function in $L^1$
- We then have
\begin{eqnarray*}
\left\Vert \left(\varphi_{\varepsilon}\ast g\right)-g\right\Vert _{1} & = & \int_{\mathbb{R}}\left|g\left(x\right)-\int_{\mathbb{R}}\varphi_{\varepsilon}\left(x-y\right)\cdot g\left(y\right)\, dy\right|\, dx\\
& \overset{\left(\ast\right)}{\leq} & \int_{\mathbb{R}}\int_{\mathbb{R}}\varphi_{\varepsilon}\left(x-y\right)\cdot\left|g\left(x\right)-g\left(y\right)\right|\, dy\, dx\\
& \overset{\text{Fubini}}{=} & \int_{\mathbb{R}}\int_{\mathbb{R}}\varphi_{\varepsilon}\left(x-y\right)\cdot\left|g\left(x\right)-g\left(y\right)\right|\, dx\, dy\\
& \overset{z=x-y}{=} & \int_{\mathbb{R}}\int_{\mathbb{R}}\varphi_{\varepsilon}\left(z\right)\cdot\left|g\left(z+y\right)-g\left(y\right)\right|\, dz\, dy\\
& \overset{\text{Fubini}}{=} & \int_{\mathbb{R}}\varphi_{\varepsilon}\left(z\right)\cdot\left\Vert g-T_{-z}g\right\Vert _{1}\, dz\\
& \overset{\left(\ast\ast\right)}{\leq} & \sup_{\left|z\right|\leq\varepsilon}\left\Vert g-T_{-z}g\right\Vert _{1}\xrightarrow[\varepsilon\downarrow0]{}0.
\end{eqnarray*}
Here, we used $\int_{\mathbb{R}}\varphi_{\varepsilon}\left(x-y\right)\, dy=1$
at the step marked with $\left(\ast\right)$ and $\int_{\mathbb{R}}\varphi_{\varepsilon}\left(z\right)\, dz=1$
as well as ${\rm supp}\left(\varphi_{\varepsilon}\right)\subset\left(-\varepsilon,\varepsilon\right)$
in the step marked with $\left(\ast\ast\right)$.
By standard results on differentiation under the integral sign, we
see that $\varphi_{\varepsilon}\ast g$ is $C^{1}$ (actually even
$C^{\infty}$) with derivative
$$
\left(\varphi_{\varepsilon}\ast g\right)'\left(x\right)=\int_{\mathbb{R}}\varphi_{\varepsilon}'\left(x-y\right)\cdot g\left(y\right)\, dy=\left(\left(\varphi_{\varepsilon}'\right)\ast g\right)\left(x\right).
$$
This is for example discussed here Differentiability of Convolutions
Now comes the interesting part, because we actually want to have $$\left(\varphi_{\varepsilon}\ast f\right)'\left(x\right)=\left(\varphi_{\varepsilon}\ast\left(f'\right)\right)\left(x\right)$$
for $f\in W^{1,1}\left(\left(0,1\right)\right)$. We will show that
this is indeed true for $x\in\left(\varepsilon,1-\varepsilon\right)$.
To see this, note that we have $\gamma_{x,\varepsilon}\in C_{c}^{\infty}\left(\left(0,1\right)\right)$
for
$$
\gamma_{x,\varepsilon}\left(y\right):=\varphi_{\varepsilon}\left(x-y\right),
$$
as $y\in{\rm supp}\left(\gamma_{x,\varepsilon}\right)$ implies $x-y\in{\rm supp}\left(\varphi_{\varepsilon}\right)\subset\left(-\varepsilon,\varepsilon\right)$
and thus $y\in B_{\varepsilon}\left(x\right)\subset\left(0,1\right)$
by our choice of $x$.
Using the definition of the weak derivative (at $\left(\ast\right)$),
we thus see
\begin{eqnarray*}
\left(\varphi_{\varepsilon}\ast f\right)'\left(x\right) & = & \int_{\mathbb{R}}\varphi_{\varepsilon}'\left(x-y\right)\cdot f\left(y\right)\, dy\\
& = & -\int_{\mathbb{R}}f\left(y\right)\cdot\gamma_{x,\varepsilon}'\left(y\right)\, dy\\
& \overset{\left(\ast\right)}{=} & \int_{\mathbb{R}}f'\left(y\right)\cdot\gamma_{x,\varepsilon}\left(y\right)\, dy\\
& = & \left(\varphi_{\varepsilon}\ast\left(f'\right)\right)\left(x\right).
\end{eqnarray*}
We know $\varphi_{\varepsilon}\ast\left(f'\right)\xrightarrow[\varepsilon\downarrow0]{L^{1}\left(\mathbb{R}\right)}f'$,
where $f'$ is extended onto all of $\mathbb{R}$ (by zero). Using
$\left(\varphi_{\varepsilon}\ast f\right)'=\varphi_{\varepsilon}\ast\left(f'\right)$
on $\left(\varepsilon,1-\varepsilon\right)$, we get
$$
\left(\varphi_{\varepsilon}\ast f\right)'\xrightarrow[\varepsilon\downarrow0]{L^{1}\left(\left(\delta,1-\delta\right)\right)}f'
$$
for any $\delta\in\left(0,\frac{1}{2}\right)$. But for $x\in\left(\delta,1-\delta\right)$,
we have
$$
F_{\varepsilon}\left(x\right)-F_{\varepsilon}\left(\frac{1}{2}\right)=\int_{1/2}^{x}F_{\varepsilon}'\left(t\right)\, dt=\int_{1/2}^{x}\left(\varphi_{\varepsilon}\ast\left(f'\right)\right)\left(t\right)\, dt\qquad\left(\dagger\right)
$$
and
$$
\left|\int_{1/2}^{x}\left(\varphi_{\varepsilon}\ast\left(f'\right)\right)\left(t\right)\, dt-\int_{1/2}^{x}\left(\varphi_{\theta}\ast\left(f'\right)\right)\left(t\right)\, dt\right|\leq\int_{\delta}^{1-\delta}\left|\left(\varphi_{\varepsilon}\ast f'\right)\left(t\right)-\left(\varphi_{\theta}\ast f'\right)\left(t\right)\right|\, dt\xrightarrow[\varepsilon,\theta\downarrow0]{}0,
$$
where we note that the expression before the limit is independent
of $x\in\left(\delta,1-\delta\right)$.
This shows that the right-hand side of $\left(\dagger\right)$ is
uniformly Cauchy, so that $F_{\varepsilon}\xrightarrow[\varepsilon\downarrow0]{}F$
uniformly on $\left(\delta,1-\delta\right)$ for some (necessarily
continuous) function $F\in C^{0}\left(\left(\delta,1-\delta\right)\right)$.
But we also know $F_{\varepsilon}\xrightarrow[\varepsilon\downarrow0]{L^{1}\left(\left(\delta,1-\delta\right)\right)}f$,
which together implies (why?) $f=F$ almost everywhere.
EDIT: As Harold pointed out in a comment, the uniform Cauchy property of $(\dagger)$ is in general not sufficient to get convergence. But in the current setting, we can fix this as follows: We know $F_\varepsilon \to f$ in $L^1$, so there is some sequence $\varepsilon_n \to 0$ with $F_{\varepsilon_n} \to f$ almost everywhere. Fix some $x_0 \in (\delta, 1-\delta)$ with $F_{\varepsilon_n}(x_0) \to f(x_0)$. We now have $$F_{\varepsilon_{n}}(x)=F_{\varepsilon_{n}}(x)-F_{\varepsilon_{n}}(1/2)-\left(F_{\varepsilon_{n}}(x_{0})-F_{\varepsilon_{n}}(1/2)\right)+F_{\varepsilon_{n}}(x_{0}),$$ where the right hand side converges uniformly by $(\dagger)$ and because $F_{\varepsilon_n}(x_0)$ converges. Now replace all limits $\varepsilon \to 0$ by "going to zero along the sequence $(\varepsilon_n)_n$".
We conclude that $F$ is a continuous version of $f$ on $\left(\delta,1-\delta\right)$
with
$$
F\left(x\right)=\lim_{\varepsilon\downarrow0}F_{\varepsilon}\left(x\right)=\lim_{\varepsilon\downarrow0}F_{\varepsilon}\left(\frac{1}{2}\right)+\int_{1/2}^{x}F_{\varepsilon}'\left(t\right)\, dt=F\left(\frac{1}{2}\right)+\int_{1/2}^{x}f'\left(t\right)\, dt.
$$
As $\delta\in\left(0,\frac{1}{2}\right)$ was arbitrary, we see
$$
f\left(x\right)=F\left(\frac{1}{2}\right)+\int_{1/2}^{x}f'\left(t\right)\, dt\qquad\left(\ddagger\right)
$$
for almost every $x\in\left(0,1\right)$.
- If $f$ is now two times weakly differentiable, the above argument
shows that $f'$ has a continuous representative, so that the right
hand side of $\left(\ddagger\right)$ is actually a $C^{1}$-function
and thus a $C^{1}$-version of $f$.
