First, solve the homogeneous equation:
$$Dy_1=y_1+2y_2$$
$$Dy_2=3y_1+4y_2$$
I am assuming you already know how to do this, so I will just write the solution here:
$$\left[\begin{matrix} y_1 \\ y_2\end{matrix}\right]=\left[\begin{matrix} \frac{1}{6}(-3-\sqrt{33})e^{(5-\sqrt{33})/2} & \frac{1}{6}(-3+\sqrt{33})e^{(5+\sqrt{33})/2} \\ e^{(5-\sqrt{33})/2} & e^{(5+\sqrt{33})/2}\end{matrix}\right]\left[\begin{matrix} c_1 \\ c_2\end{matrix}\right] \text{ where } c_1,c_2\in\mathbb{R}$$
Since there is one negative eigenvalue $\frac{5-\sqrt{33}}{2}$ and one positive eigenvalue $\frac{5+\sqrt{33}}{2}$, the origin is a saddle point and the system is unstable. Therefore, at the very least, you now know the stability of the ODE system.
Now, to actually solve the non-homogeneous equation, use variation of parameters by changing $c_1,c_2$ to $u_1(t),u_2(t)$:
$$\left[\begin{matrix} y_1 \\ y_2\end{matrix}\right]=\left[\begin{matrix} \frac{1}{6}(-3-\sqrt{33})e^{(5-\sqrt{33})/2} & \frac{1}{6}(-3+\sqrt{33})e^{(5+\sqrt{33})/2} \\ e^{(5-\sqrt{33})/2} & e^{(5+\sqrt{33})/2}\end{matrix}\right]\left[\begin{matrix} u_1(t) \\ u_2(t)\end{matrix}\right]$$
$$\left[\begin{matrix} \frac{1}{6}(-3-\sqrt{33})e^{(5-\sqrt{33})/2} & \frac{1}{6}(-3+\sqrt{33})e^{(5+\sqrt{33})/2} \\ e^{(5-\sqrt{33})/2} & e^{(5+\sqrt{33})/2}\end{matrix}\right]\left[\begin{matrix} u_1'(t) \\ u_2'(t)\end{matrix}\right]=\left[\begin{matrix} \frac{3}{x^4} \\ \frac{3}{x^4}\end{matrix}\right]$$
Here, you can use the second equation to solve for $u_1'(t),u_2'(t)$, integrate to solve for $u_1(t),u_2(t)$, and then plug into the first equation to solve for $y_1,y_2$. Good luck!
To remain second order in the space approximation, you would use the central difference quotient. You would then subtract a linear approximation of the right side from both sides, so that you get an iterative fixed-point scheme $A\vec y^{\rm new}=G(\vec y^{\rm old})$, where $\vec y$ is the vector of sample points. Ideally this would be Newton-like so that you get super-linear convergence.
If you want to generalize this scheme you get to a multiple shooting approach with a collocation method, which is for first order systems, so set $v=y'$, $v'=f(v,y,x)$. Staying in second order this would be the implicit trapezoidal method,
\begin{align}
\frac{y_{i+1}-y_i}{\Delta x}&=\frac{v_{i+1}+v_i}2\\
\frac{v_{i+1}-v_i}{\Delta x}&=\frac{f(v_{i+1},x_{i+1},v_{i+1})+f(v_i,y_i,x_i)}2\\
\end{align}
You could also try to base your system on Numerov's method (wiki)
Best Answer
I'm not sure what you mean by "applying this method to the system... You have written down the linear system to be solved (apart from the correction mentioned by @Ian), now you just need to solve it. Your equations are $$ -3y_0 + 4y_1-y_2 = - 6h $$ $$ y_{i-1}-(2+h^2)y_i+y_{i+1} = 8h^2+3h^2 x_i(1-x_i), \quad i = 1, \cdots, n $$ $$ y_{n-1}-4y_n+3y_{n+1} = 2h (e-e^{-1}+3) $$
so, apart from the first and last lines, you have a tridiagonal system. Below you can visualize the numerical solution with $n=10$ (dots), compared to the exact solution $y(x)=-2 + e^{-t} + e^t - 3 t + 3 t^2$ (line).
It is also very instructive to observe the behaviour of the numerical solution when you use first order approximations of the derivatives to enforce the Neumann conditions, i.e. what happens if you replace the fisrt and last equations by $y_1-y_0=-3h$ and $y_{n+1}-y_n = h(e-e^{-1}+3)$.