Finding the general solution to $-u''(x)=f(x)$ for general $f(x)$ 
Problem: Use variation of parameters to verify that $$u(x)=\frac{L-x}{L}\left(\gamma_0+\int_{0}^{x} yf(y) \ dy\right)+\frac{x}{L}\left(\gamma_L+\int_{x}^{L}(L-y)f(y) \ dy\right)  \ \ \text{for} \ 0\leq x\leq L$$ solves $-u''(x)=f(x)$ for $0<x<L$ with $u(0)=\gamma_o$ and $u(L)=\gamma_L$ for a general  $f(x)$. Let $$u_1(x)=L-x \ \ \text{and} \ \ u_2(x)=x$$ and write $$u(x)=v_1(x)u_1(x)+v_2(x)u_2(x).$$

I proceed via variation of parameters, where
$$v_1'(x)=-\frac{u_2(x)\left(-f(x)\right)}{W(x)}=\frac{u_2(x)f(x)}{W(x)} \ \ \text{and} \ \ v_2'(x)=\frac{u_1(x)\left(-f(x)\right)}{W(x)}=-\frac{u_1(x)f(x)}{W(x)},$$ noting that $u''(x)=-f(x)$. Here, $W(x)$ denotes the Wronskian such that $W(x)=L$. Thus $$v_1'(x)=\frac{xf(x)}{L} \ \ \text{and} \ \ v_2'(x)=\frac{(L-x)f(x)}{L}$$ such that $u(x)=v_1 (x)\left(L-x\right)+v_2(x)x$. To find $v_1(x)$ and $v_2(x)$ I integrate:
$$v_1(x)=\frac{1}{L}\int_{0}^{x} yf(y) \ dy\ +C_1 \ \text{and} \ \ v_2(x)=\frac{1}{L}\int_{0}^{x}(L-y)f(y) \ dy +C_2 \ \ \ C_1,C_2\in\mathbb{R}.$$ I do not understand why the bounds of integration for $v_2(x)$ are between $x$ and $L$. After applying the boundary conditions, I find that $$u(x)=\frac{L-x}{L}\left(\gamma_0+\int_{0}^{x} yf(y) \ dy\right)+\frac{x}{L}\left(\gamma_L+\int_{0}^{x}(L-y)f(y) \ dy\right),$$ which does not match what is provided.
 A: There seems to be some issues with your parameters : $v_1'(x)$ and $v_2'(x)$.
I would try to give a brief procedure how to go about solving the problem:


*

*Find the homogeneous solution to the equation $-u''(x) = f(x)$: 


$$u''(x) = 0$$
$$\implies u(x) = A + Bx$$
A and B are constants. The homogeneous solutions are : $u_1(x) = 1$ and $u_2(x) = x$.


*

*Find the parameters:
$$W(u_1, u_2) = u_1u_2' - u_1'u_2 = 1\cdot 1 - 0\cdot x = 0$$
$$v_1'(x) = \frac{u_2(x)f(x)}{W(u_1, u_2)}; \hspace{10mm}  v_2'(x) = \frac{-u_1(x)f(x)}{W(u_1, u_2)}$$
$$\implies v_1'(x) = {xf(x)}; \hspace{10mm} v_2'(x) = -f(x)$$
$$\text{Now, } v_1(x) = \int_0^x yf(y) dy + C_1; \hspace{10mm} v_2(x) = -\int_0^x f(y) dy + C_2$$
where, $c_1$ and $C_2$ are constants.


*

*Solve for the particular solution: 
$$u_p(x) = v_1(x) u_1(x) + v_2(x) u_2(x) = v_1(x) + v_2(x) x$$

*The general solution is:
\begin{align}
   u(x) &= A + Bx + u_p(x) \\
         &= A + Bx + \int_0^x yf(y) dy + C_1 + x\left(-\int_0^x f(y) dy + C_2\right) \\
          &= (A+C_1) + (B+ C_2) x + \int_0^x yf(y) dy - x\int_0^x f(y) dy \\
          &= D + Ex + \int_0^x yf(y) dy - x\int_0^x f(y) dy
  \end{align}

*Apply Boundary conditions:
$$u(0) = \gamma_0,u(L) = \gamma_L$$
\begin{align}
     u(x = 0) = \gamma_0 &\implies D = \gamma_0 \\
     u(x = L) = \gamma_L &\implies \gamma_0 + EL + \int_0^L yf(y) dy - L\int_0^L f(y) dy = \gamma_L \\
   &\implies E = \frac{\gamma_L}{L} - \frac{\gamma_0}{L} - \frac{1}{L}\int_0^L yf(y) dy + \int_0^L f(y) dy
 \end{align}

*Subsitute the expressions for $D$ and $E$ to get particular solution:
(and for terms with no coefficients and/or with only $x$, multiply and divide by $L$).
\begin{align}
 u(x) &= \gamma_0 + x\left[\frac{\gamma_L}{L} - \frac{\gamma_0}{L} - \frac{1}{L}\int_0^L yf(y) dy + \int_0^L f(y) dy\right] + \int_0^x yf(y) dy - x\int_0^x f(y) dy \\
      &= \frac{L}{L}\gamma_0 + \frac{x\gamma_L}{L} - \frac{x\gamma_0}{L} - \frac{x}{L}\int_0^L yf(y) dy + \frac{xL}{L}\int_0^L f(y) dy + \frac{L}{L}\int_0^x yf(y) dy - \frac{xL}{L}\int_0^x f(y) dy \\
     &= \frac{L-x}{L}\gamma_0 + \frac{x\gamma_L}{L} - \frac{x}{L}\int_0^L yf(y) dy  + \frac{xL}{L}\left[\int_0^L f(y) dy - \int_0^x f(y) dy\right] + \frac{L}{L}\int_0^x yf(y) dy\\
  &= \frac{L-x}{L}\gamma_0 + \frac{x\gamma_L}{L} - \frac{x}{L}\left[\int_0^x yf(y) dy + \int_x^L yf(y) dy\right]  + \frac{xL}{L}\int_x^L f(y) dy + \frac{L}{L}\int_0^x yf(y) dy \\
&= \frac{L-x}{L}\gamma_0 + \frac{L-x}{L}\int_0^x yf(y) dy + \frac{x}{L}\gamma_L + \frac{x}{L}\int_x^L (L-y)f(y) dy \\
\implies  u(x)  &= \frac{L-x}{L}\left(\gamma_0 + \int_0^x yf(y) dy\right) + \frac{x}{L}\left(\gamma_L + \int_x^L (L-y)f(y) dy \right)
\end{align}
Hope it helps!
A: You just made a small mistake in calculating $C_2$. Since
\begin{eqnarray}
u(x)&=&u_1(x)v_1(x)+u_2(x)v_2(x)\\
&=&\frac{L-x}{L}\bigg[\int_{0}^{x} yf(y) \ dy\ +C_1\bigg]-\frac{x}{L}\bigg[\int_{0}^{x}(L-y)f(y) \ dy +C_2\bigg]
\end{eqnarray}
using $u(0)=\gamma_0,u(L)=\gamma_L$, one has
$$ C_1=\gamma_0$$
and
$$ \int_{0}^{L}(L-y)f(y) \ dy +C_2=\gamma_L. $$
The latter gives
$$ C_2=-\int_{0}^{L}(L-y)f(y) \ dy+\gamma_L. $$
So
\begin{eqnarray}
u(x)&=&\frac{L-x}{L}\bigg[\int_{0}^{x} yf(y) \ dy\ +\gamma_0\bigg]-\frac{x}{L}\bigg[\int_{0}^{x}(L-y)f(y) \ dy -\int_{0}^{L}(L-y)f(y) \ dy+\gamma_L\bigg]\\
&=&\frac{L-x}{L}\bigg[\int_{0}^{x} yf(y) \ dy\ +\gamma_0\bigg]+\frac{x}{L}\bigg[\int_{x}^{L}(L-y)f(y) \ dy+\gamma_L\bigg].
\end{eqnarray}
