What is the process for deducing Galerkin approximation for an elliptic boundary value problem? For example, for a problem like--
$$ -\Delta u + cu = f\ \ \mathrm{in}\ Ω $$
$$ u = g\ \ \mathrm{on}\ \ \Gamma_1 $$
$$ \frac{∂u}{∂n} + pu = h\ \ \mathrm{on}\ \ \Gamma_2 $$
here, $∂Ω = \Gamma_1 \cup \Gamma_2$ and $\Gamma_1\cap \Gamma_2 = \phi$ and $c \geq 0.$
$f, g, h, p$ are sufficiently smooth functions.
I’ve calculated the weak solution of this problem to be
$$ \int_\Omega (\nabla u\cdot \nabla \phi + cu\phi)dw + \int_{\Gamma_2} (pu-h)\phi ds = \int_{\Omega} (f\phi) dw $$
here $\Omega, \Gamma_2$ with the integral sign $\int$ are the domain identifier.
So how can I propose a galerkin approximation from this weak form of the PDE
The idea of the Galerkin FEM is to approximate this solution. After the steps you have mentioned, just take a finite element subspace $V_h\subset V=H^1_0{\Omega}$. Then, as $V_h$ is finite, it will have a finite basis $\lbrace \phi\rbrace_{i=1}^n$ where $n$ is the dimension of the space. Then you can approximate the solution $u$ by $u_h$ where $u_h$ is a linear combination of $\lbrace \phi\rbrace_{i=1}^n$, i.e, $u_h=\sum_{i=1}^n\phi_iu_i$. Now using, $\phi=\phi_i$ in your equations, you will get a system of equations of the form $Au=b$ where $A$ is referred to as stiffness matrix and $b$ as load vector. You can refer to any standard FEM book for more details.
There is one small thing, you might be confused because of non-homogenous BCs. In this case the Dirichlet BCs are incorporated in the function space and the Neumann conditions will come in the weak form (as you can see in your defintion).