How to get the linear equation system for finite element method from the variational formulation Let the problem be
$$-u'' + a(x) u = f , \;x \in \Omega = ]0,1[ , \;u(0) = \alpha ; \;u(1) = \beta,$$
where $f \in L^2(\Omega) , a(x) \geq a_0 > 0 , a(x) \in L^{\infty}(\Omega).$
This problem admits a unique solution in $V = H^1(\Omega).$
My question is: 

How we can use the finite elements $\mathbb{P}_1$ to prove that this problem, in the finite element space, is equivalent to a linear system $A U = b$ .

BTW: for this, we can use the relation for linear function $\psi$:
 $$\dfrac{1}{h} \int_{x_i}^{x_{i+1}} \psi(x) dx = \psi\left(\dfrac{x_i + x_{i+1}}{2}\right).$$
 A: The nodal basis function for $\mathbb{P}_1$-elements is exactly as you said in the comments, a.k.a. hat functions (see Figure 11.2 here):
$$
w_j(x_i) = \begin{cases}1\quad \text{ when }i = j,
\\
0\quad \text{ otherwise },\end{cases}
$$
and the support of a hat function $w_j$ is $(x_{j-1},x_{j+1})$ in which the hat function is linear. First we assume the finite element solution 
$$u_h(x) = \alpha w_0(x) + \beta w_n(x) + \sum^{n-1}_{n=1} u_i w_i(x),\tag{1}$$
so that $u_h$ satisfies the boundary condition, and $u_i$'s above are coefficients (degrees of freedom sometimes we call them) we want to find. Now we wanna plugging the hat function one by one in the weak formulation:
$$
\int_0^1 \left(u'(x)v'(x) + a(x)u(x)v(x)\right)\,dx = \int_0^1 f(x) v(x)\,dx.\tag{2}
$$
Plugging (1) into (2), then replacing $v(x)$ by hat functions $w_j$, $j = 1,\ldots,n-1$ (do not include the hat functions on the end points for the test function space has zero boundary condition), we have:
$$
\int_0^1 \left(\big(\sum^{n-1}_{n=1} u_i w_i(x)\big)'w_j'(x) + a(x)\big(\sum^{n-1}_{n=1} u_i w_i(x)\big)w_j(x)\right)\,dx = \int_0^1 f(x) w_j(x)\,dx,$$
pulling the summation outside the integral for $u_i$'s are numbers:
$$
\sum^{n-1}_{n=1} u_i \int_0^1 \left(w'_i(x)w_j'(x) + a(x) w_i(x) w_j(x)\right)\,dx = \int_0^1 f(x) w_j(x)\,dx,
$$
Now we have a linear system:
\begin{gather}
&AU = b, \quad \text{ where }
\\
&A_{ji} =  \int_0^1 \left(w'_i(x)w_j'(x) + a(x) w_i(x) w_j(x)\right)\,dx : \text{ an } (n-1)\times (n-1) \text{ matrix } ,
\\
&U_{i} = u_i:  \text{ an } (n-1)\text{-column vector } ,
\\
&b_{j} = \int_0^1 f(x) w_j(x)\,dx: \text{ an } (n-1)\text{-column vector } .
\end{gather}
Now the last task is to evaluate:
$$
A_{ji} =  \int_0^1 \left(w'_i(x)w_j'(x) + a(x) w_i(x) w_j(x)\right)\,dx,
$$
and $A_{ji}$ is zero if two hat functions do not overlap. $ w_i(x)$ and $ w_j(x)$ overlap when $i= j,j-1,j+1$, so the rest work is evaluate:
$$
\begin{aligned}
&A_{ii} =  \int_{x_{i-1}}^{x_{i+1}} \left(w'_i(x)w_i'(x) + a(x) w_i(x) w_i(x)\right)\,dx,
\\
&A_{i,i-1} =  \int_{x_{i-2}}^{x_{i+1}} \left(w'_{i-1}(x)w_i'(x) + a(x) w_{i-1}(x) w_i(x)\right)\,dx,
\\
&A_{i+1,i} =  \int_{x_{i-1}}^{x_{i+2}} \left(w'_{i}(x)w_{i+1}'(x) + a(x) w_{i}(x) w_{i+1}(x)\right)\,dx.
\end{aligned}\tag{3}
$$
I suggest you try this yourself with the integration formula you gave for linear functions.

UPDATE 1: Nodal basis "hat" functions are:
$$
w_i(x)=\frac{x−x_{i−1}}{h} \;\text{ in } [x_{i−1},x_i], \quad \text{and}\quad  w_i(x)=−\frac{x−x_{i+1}}{h} \;\text{ in } [x_i,x_{i+1}].
$$
First integration in $A_{ii}$, notice here we assume all $x_i$'s are evenly spaced, $x_i - x_{i-1} = x_{i+1}-x_i = h$:
$$
\int_{x_{i-1}}^{x_{i+1}} w'_i(x)^2 = \frac{x_i−x_{i−1}}{h^2} + \frac{x_{i+1}−x_{i}}{h^2} = \frac{2}{h}.
$$
For the second integration, we can't do much, unless we assume $a$ is a piecewise constant, otherwise we have to use quadrature rule, and you can't pull $w_i^2$ out. Assume $a= a_{i-1}$ in $[x_{i−1},x_i]$, and $a = i+1$ in $[x_i,x_{i+1}]$
$$
\int_{x_{i-1}}^{x_{i+1}}a(x) w_i(x)^2 = \frac{a_{i-1}h}{3} + \frac{a_i h}{3} = \frac{h(a_{i-1}+a_i)}{3}.
$$
For $A_{i,i-1}$ and $A_{i,i+1}$, beware of the overlapping part of different hat functions, so $A_{i,i-1}$ is really the integration from $x_{i-1}$ to $x_i$, and $A_{i,i+1}$ is really the integration from $x_{i}$ to $x_{i+1}$ .

UPDATE 2: The numerical integration of an arbitrary polynomial on interval $[\alpha,\beta]$ can use the following handy rules for a single term:
$$
\int^{\alpha}_{\beta} w_{\alpha}(x)^m w_{\beta}(x)^n\,dx = \frac{m!\,n!}{(m+n+1)!} |\beta-\alpha|,
$$
where $w_a(x)$  is 1 at $a$ and 0 at $b$ and linear within, similarly $w_b$ is the hat function for $b$, in your case, if coefficient is constant:
$$
\int_{x_{i-1}}^{x_{i}}  w_i(x)^2\,dx = \int_{x_{i-1}}^{x_{i}}  w_i(x)^2 w_{i-1}(x)^0\,dx= \frac{2!\,0!}{(2+0+1)!} |x_{i} -x_{i-1} | = \frac{1}{3}h,
$$
and
$$\int_{x_{i-1}}^{x_{i}}  w_i(x) w_{i-1}(x) \,dx= \frac{1!\,1!}{(1+1+1)!} |x_{i} -x_{i-1} | = \frac{1}{6}h.$$
