Question: I was reading about physics and came across the method of using separation of variables to solve specific PDEs, but I can't figure out why the specific solutions give rise to the general solution (the book didn't give any explanation for all these).

The specific example in the book was the Laplace Equation in $2$ variables:
$$\frac {\partial^2 V}{\partial x^2}+\frac {\partial^2 V}{\partial y^2}=0$$
For the above example, separation of variable is essentially solving for the eigen-vectors of the operator $\frac {\partial^2 }{\partial x^2}$ and $\frac {\partial^2 }{\partial y^2}$, which are Hermitian and commutes with each other. I know that in the finite dimensional case, such operators are simultaneously diagonalizable, then solving for the eigen-vectors will give all the solution, but I'm not sure does this work for infinite dimension. I'm also not sure does this approach works in the general case, for other PDEs that can be solved by separation of variable.

All the other post I find on here are all explaining how or when separation of variable work, instead of why such techniques will give the general solutions.

Another side question is: What kind of classes will cover these topics? The only undergraduate class that seems relevant at my university is Linear Analysis, which doesn't cover this. The graduate PDE sequence have graduate Real Analysis sequence as pre-requisite, which I don't think I'll be able to take soon.

Answer: 
There are several key ingredients I will briefly describe here. I won't go into too much detail as you've mentioned that you don't have a graduate real analysis background yet. But indeed a full description of the theory is a standard part of a graduate course in linear PDE. So I hope that answers your side question as well.

1. We start with a strongly elliptic linear operator (such as the Laplacian) and, along with some nice boundary condition, we restrict to some appropriate solution (Hilbert) space.
2. In that solution space, we can prove under fairly general conditions that the eigenvalues of the operator are countable and that eigenvectors (eigenfunctions) form an orthogonal basis for the solution space. This is the infinite-dimensional generalization of the diagonalizability result from regular matrix theory. The proof relies on the spectral theorem for compact operators. The key here is that, up to a shift, the inverse of a strongly elliptic operator is compact.
3. This demonstrates that if we can construct all the eigenvectors of the operator, the general solution can be written as a decomposition of these eigenvectors.
4. It remains to find the eigenvectors; in special cases (most famously, 2D Laplacian on a rectangle) this can be done via separation of variables. Therefore it remains to address "Why does separation of variables produce all eigenvectors?" To answer this question, we note that we proved that the eigenvectors form a complete basis. Next, we see that because of the specific symmetry of the Laplacian on the rectangle, using separation of variables reduces the problem to a pair of second-order equations in one-dimension; in this process we produce the eigenvectors of these one-dimensional operators, and then from the existing theory (in particular, Sturm-Liouville theory) we know that we have produced a set of functions that span the space. As we have produced a basis, no other eigenvectors are needed to form a general solution.
The answer by @Christopher is very complete and definitely better than what this answer will be. But I would like to make some comments on Separation of Variables.

Separation of Variables is a process of splitting a multi-dimensional problem, into several single dimensional problems. However, this relies on an inherent symmetry of the domain, which itself determines the coordinates which allow for separation of variables.

If the question is posed in a rectangle, then it is quite natural that the problem given in rectangular coordinates can be broken down into two one dimensional problems in each orthogonal dimension. If the problem is posed on a circle, then polar coordinates are required. However, if the problem is given on a completely arbitrary domain then it is unlikely that you could find a coordinate system that can reflect the symmetry of the domain and allow for separation of variables.

If you get deeper into Lie theory, one can describe a group theoretic method of determining the possible coordinate systems that allow for a given equation to be separable. However, I don't think I have a deep enough understanding on this to comment further.