Numerical solution of $y′′+2y=−x$? How to solve $y′′+2y=−x$ differential equation numerically. 
$y′(1)=0$ and $y(0)=0$ ?
 A: For example, you can apply a 2nd order finite difference scheme to your ODE with the following approximations over the grid $x_i = x_0 + i \Delta x$, $i = 0, \ldots, N$, $\Delta x = (x_N-x_0)/\Delta x$, $x_0 = 0, x_N = 1$:
$$y'' \approx \frac{Y_{i-1} - 2 Y_i + Y_{i+1}}{\Delta x^2},$$ where $Y_i \approx y(x_i)$. Plug this into the original ODE to have:
$$Y_{i-1} +2 (\Delta x^2 -1) Y_i + Y_{i+1} = - x_i \Delta x^2.$$
Let's not regard the limits of $i$ for which this equation is valid. Let's discretize the boundary conditions first. Thus:
$$y(0)=0\Rightarrow Y_0 = 0, \quad y'(1) = 0 \Rightarrow \frac{Y_{N+1}-Y_{N-1}}{2\Delta x} = 0,$$ so the last equation introduce a fictitious node. This allows us to write that: $Y_{N+1} = Y_{N-1}$, so if we introduce this information into the finite difference equation for $i =1$ and $i =N$ (which result to be the limits of validity) we come up to:


*

*$i = 1$: The equation becomes: $2 (\Delta x^2 -1) Y_1 + Y_{2} = - x_i \Delta x^2.$

*$i = N$: The equation becomes: $2Y_{N-1} +2 (\Delta x^2 -1) Y_N = - x_i \Delta x^2.$


And with this information we can conclude that the problem given by the finite difference equations, with its resepective modifications, can be arranged in matrix form as follows:
$$\left(\! \!
\begin{array}{ccccc}
 2\left(\text{$\Delta $x}^2-1\right) & 1 & \cdots  & 0 & 0 \\
 1 & 2\left(\text{$\Delta $x}^2-1\right) & \cdots  & 0 & 0 \\
 \vdots  & \vdots  & \ddots & \vdots  & \vdots  \\
 0 & 0 & \cdots  & 2\left(\text{$\Delta $x}^2-1\right) & 1 \\
 0 & 0 & \cdots  & 2 & 2\left(\text{$\Delta $x}^2-1\right)
\end{array}\! \!
\right)\left(\! \!
\begin{array}{c}
 Y_1 \\
 Y_2 \\
 \vdots  \\
 Y_{N-1} \\
 Y_N
\end{array}\! \!
\right)=\text{$\Delta $x}^2\left(\! \!
\begin{array}{c}
 x_1 \\
 x_2 \\
 \vdots  \\
 x_{N-1} \\
 x_N
\end{array}\! \!
\right),$$ and solve for the vector of unknowns, $Y_i$.
Cheers!

Edit: Matlab code for solving this problem (feel free to change the number of intervals, $N$):
N = 100;
dx = 1/N;
x = 0:dx:1;

M = zeros(N,N);

diagonal = 2*(dx^2-1)*ones(N,1);
sup_diag = ones(N-1,1);
inf_diag = ones(N-1,1);
inf_diag(N-1) = 2;

M  = diag(diagonal,0) + diag(sup_diag,1) + diag(inf_diag,-1);

C = zeros(N,1);
C = dx^2*x(2:N+1)';

Y = zeros(N+1,1);
Y(2:N+1) = M\C;

plot(x,Y);

A: Lots of numerical methods can solve this equation.
For example, FDM is a simple way to solve it numerically. The numerical scheme is established as follow:
Discretizing the $[0,1]$ by $\{0=x_0,x_1,x_2,\dots,x_N=1\}$ where $x_i$ ($i=0\dots N$) is the set of partition points and $N$ is partition steps. Let $\{y_0,y_1,\dots,y_N\}$ be the finite different approximation of $y(x)$ at $\{x_i\}$. We use the middle point formula to approximate $y''$ and then obtain the following numerical scheme of original equation
\begin{equation}
\frac{y_{i+1}-2y_i+y_{i-1}}{h_i^2}+2y_i=-x_i\quad i=1\dots N-1
\end{equation}
where $h_i=x_{i+1}-x_i$ is the length of partition interval. Note that only $y_i$ in these equations are unknown. This scheme gives $N-1$ equations but $N+1$ unknown variables. Note $y(0)=0$ and $y'(1)=0$ also give two conditions:
\begin{equation}
y_0=0\\
\frac{y_N-y_{N-1}}{h_{N-1}}=0
\end{equation}
Adding these two equations into previous systems, we obtain a linear systems have $N+1$ variables and $N+1$ equations. Also, we can show the system is unique solvable. So the numerical result is given by the solution of the linear system.
What is more, other numerical methods, such as FEM and FVM etc. can also solve this problem. You can find them in references about numerical method of PDE and ODE.
