Derivation of the perihelion advance relation
Let's start with this equation
$$\frac{d^2 u}{d\phi^2}+u = \frac{M}{L^2}+3Mu^2, \tag 1$$
and let's go step-by-step. Remember that this equation is valid in units with $G=c=1$.
Step 1. First of all, since this is a perturbative analysis, we have to consider a circular orbit with a constant radius ($r_0$). This ($u_0=1/r_0$) obviously satisfies Eq. (1), i.e.,
$$u_c = \frac{M}{L^2}+3Mu_0^2, \tag 2$$
in which $\frac{d^2 u_0}{d\phi^2}=0$. You have made a mistake here, in fact, these derivatives are zero: $\frac{d^2 u_0}{d\phi^2}=\frac{d u_0}{d\phi}=0$.
Step 2. Next, we have to consider a slightly non-circular orbit as
$$u(\phi, \epsilon)=u_0+\epsilon u_0 u_1=u_0(1+\epsilon u_1), \tag 3$$
which is a perturbation of the previous circular orbit in Step 1. By inserting Eq. (3) into Eq. (1) together with the Eq. (2), you will obtain the following equation
$$\epsilon \frac{{{d^2}{u_1}}}{{d{\phi ^2}}} + \epsilon {u_1} = 6\epsilon M{u_0}{u_1} + 3{\epsilon ^2}M{u_0}u_1^2, \tag {4-A}$$
or equivalently
$$\frac{{{d^2}{u_1}}}{{d{\phi ^2}}} + {u_1} = 6 M{u_0}{u_1} + 3{\epsilon}M{u_0}u_1^2. \tag {4-B}$$
Now, since it's been assumed that $\epsilon u_1 \ll 1$, the second term on the R.H.S. of the Eq. (4) is much smaller than the first term, so ignore it. This results in
$$\frac{{{d^2}{u_1}}}{{d{\phi ^2}}} + \left( {1 - 6M{u_0}} \right){u_1} = 0. \tag 5$$
Step 3. Finally, we have ended up with the harmonic oscillator equation, Eq. (5). It's solution is simple, i.e.,
$${u_1} = A\cos \left( {\omega \phi + {\phi _0}} \right), \tag 6$$
where
$$\omega = \sqrt {1 - 6M{u_0}}. \tag 7$$
I think that you know the rest of typical small steps, in summary: a) when $u_1$ is a minimum, the orbit’s perihelion occurs b) then consider a change of the argument by $2\pi$ c) next, using the binomial approximation and assuming that $u_1 M \ll 1$, yielding
$$\omega \Delta \phi = 2\pi \Rightarrow \Delta \phi = 2\pi + 6\pi M{u_0}.$$
This is the famous relation of perihelion advance in general relativity.
Applying to the planets' orbits in the solar system
This derivation was performed in units with $G=c=1$. It is easy to convert this into SI units by performing this substitution: $M \to \frac{GM}{c^2}$. So, using this substitution together with $u_0=1/r_0$, one finds
$${\rm{perihelion}} \, {\rm{advance}}\,{\rm{relation:}}\,6\pi M{u_0} \to \frac{{6\pi GM}}{{{r_0}{c^2}}},$$
in which $r_0$ is the mean orbital radius (the mean distance from Sun). Considering Sun with the mass $1.989 \times 10^{30}$ kg and the Mercury with the mean distance from Sun about $5.79 \times 10^7$ km, one finds the famous result of $43$ $\frac{\rm{arc-seconds}}{\rm{century}}$. (data from Solar System Data)
$=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!$
References
[1] T.A. Moore, A General Relativity Workbook, University Science Books (2012).
[2] S.M. Carroll, Spacetime and geometry: an introduction to general relativity, Addison-Wesley, San Francisco, USA (2004).
[3] M.P. Hobson, G.P. Efstathiou, and A.N. Lasenby, General relativity: an introduction for physicists, Cambridge University Press (2006).