Basic Definitions:
If you don't know the definitions of the basic computational complexity classes well, I strongly recommend going through Watrous' lecture. We won't be using the quantum Turing machine formulation here, unlike the formal rigorous proof by Berstein and Vazirani. Anyway, I'm including a brief discussion on the definitions here.
Polynomial-time quantum algorithms: In the quantum circuit model, a quantum algorithm $Q$ is described by a family of quantum circuits $$Q = \{Q_n:n\in \Bbb N\}.$$ Such a family is polynomial-time uniform if there exists a classical algorithm that produces a description of $Q_n$ for each input $n\in \Bbb N$, in time polynomial in $n$.
To run this algorithm on an input $x$ of length $n$, we apply $Q_n$ to $|x\rangle$, and measure the output in the standard basis:
Shorthand: Write $Q(x)$ to denote the outcome, which is a random variable in general.
Promise Problem : It is a pair $$\mathcal{A}=\{\mathcal{A}_{\text{yes}},\mathcal{A}_{\text{no}}\}$$
where $\mathcal{A}_{\text{yes}}$ and $\mathcal{A}_{\text{no}}$ are
disjoint sets of inputs (yes inputs and no inputs). Strings
not contained in $\mathcal{A}_{\text{yes}} \cup
\mathcal{A}_{\text{no}}$ are "don't care" inputs.
$\mathsf{BQP}$: The class of promise problems $\mathcal{A} = (\mathcal{A}_{\text{yes}}, \mathcal{A}_{\text{no}})$ for which there exists a polynomial time-uniform family
$$\mathrm{Q} = \{\mathrm{Q}_n: n\in \Bbb N\}$$
of quantum circuits such that, for all input strings $x$, we have:
$$x\in\mathcal{A}_{\text{yes}} \implies \text{Pr}[\mathrm{Q}(x)=1]\geq \frac{2}{3}$$
$$x\in\mathcal{A}_{\text{no}} \implies \text{Pr}[\mathrm{Q}(x)=0]\geq \frac{2}{3}$$
$\mathsf{PSPACE}$: The class of promise problems solvable in polynomial space on a deterministic Turing machine.
Source: John Watrous - Quantum Complexity Theory (Part 1) - CSSQI 2012
Note: In the gate model, a quantum algorithm is essentially a set of unitary gate operations followed by measurement. Think of $Q(x)$ as a sequence of gate operations performed on an $n$-qubit state $|x\rangle$, with measurement performed on one specific qubit at the end. As a qubit can only collapse to the state $|0\rangle$ or $|1\rangle$ when measured in the standard basis, you'll get either $0$ or $1$ as the output i.e. either $Q(x) = 1$ or $Q(x) = 0$.
Explanation:
Why is it that $\mathsf{BQP}\subseteq \mathsf{PSPACE}$? Here's an intuitive outline of the proof (a rigorous proof is left to the references in 'History and further reading'). Suppose we have an $n$ qubit quantum computer, and do a computation involving a set of $p(n)$ gates, where $p(n)$ is some polynomial in $n$. Supposing the quantum circuit starts in the state $|0\rangle$ we will explain how to evaluate in polynomial space on a classical computer the probability that it ends in the state $|y\rangle$. Suppose the gates that are executed on the quantum computer, are in the order $U_1, U_2, ..., U_{p(n)}$. Then the probability of ending up in the state $|y\rangle$ is the modulus squared of $$\langle y|U_{p(n)} ... U_2U_1|0\rangle. \tag{4.86}$$
This quantity may be estimated in polynomial space on a classical computer. The basic idea is to insert the completeness relation $\sum_x|x\rangle\langle x|=I$ between each term in $(4.86)$ obtaining
$$\langle y| U_{p(n)}...U_2U_1 |0\rangle = \sum_{x_1,...,x_{p(n)-1}} \langle y |U_{p(n)}|x_{p(n)-1}\rangle \langle x_{p(n)-1}|U_{p_{n-1}}... U_2|x_1\rangle \langle x_1|U_1|0\rangle.$$
Nielsen and Chuang prove that if your quantum circuit begins with all qubits initialized in the state $|0\dots\rangle$ (all qubits in state $|0\rangle$; N&C has labeled that state as $|0\rangle$ so don't get confused), you can evaluate in polynomial space (i.e. the memory space required on a classical computer can be written down as a polynomial in $n$) the probability that the final would be any arbitrary state $|y\rangle$.
In other words, given a certain composite quantum state $|x\rangle$, you can calculate in polynomial space the probability of a specific qubit returning $|0\rangle$ (or $|1\rangle$), upon the final measurement. You may be thinking that even the generating the description of $Q_n$'s, the $Q_n$'s themselves, and storing the results of their operations on the composite quantum state would occupy space. Yes, that is true. However, that too will only be a polynomial amount of space, as in polynomial time only polynomial space can be accessed. C.f. Polynomial-time algorithm with exponential space is eligible?
In the background, it's really similar to the path integral approach (you may skip this and jump to the exercise at the end of answer, in the first reading).
In terms of computational complexity, the Schrodinger picture
$\alpha_x|x\rangle$ and Heisenberg's density matrix $\rho$, both lead
to exponential-space simulation since we need to calculate whole
evolution of state vectors. On the other hand, the Feynman's path
integral, summing up all the histories, leads to a polynomial-space
simulation. By writing each final amplitude as a sum of contributions
from all possible paths, we can calculate the sum in
$\mathsf{PSPACE}$.
For example, the calculation of $H\otimes H|0\rangle$ can be viewed as
follows in Feynman's path integral. We calculate amplitude for each
path separately which needs polynomial space only.
In short, the point is that you can evaluate the result of your decision problem in polynomial space by summing over the desired paths instead of having to carry out those exponential size matrix multiplications. For arbitrary inputs $|x\rangle$, you'd have to decompose it into the individual standard basis vectors and carry out the calculation for each $2^n$ basis vector separately (erasing previous data as you proceed). Note the emphasized part in the following paragraph (from N&C) carefully, in case you're worried that there could be exponentially many paths to sum over.
Given that the individual unitary gates appearing in this sum are operations such as the Hadamard gate, CNOT, and so on, it is clear that each term in the sum can be calculated to high accuracy using only polynomial space on a classical computer, and thus the sum as a whole can be calculated using polynomial space, since individual terms in the sum can be erased after being added to the running total. Of course, this algorithm is rather slow, since there are exponentially many terms in the sum which need to be calculated and added to the total; however, only polynomially much space is consumed, and thus $\mathsf{BQP} \subseteq \mathsf{PSPACE}$, as we set out to show.
Exercise:
Calculate the probability of obtaining $|0\rangle$ (upon measurement), in the first (top) qubit, for this circuit using the summation over paths methods explained above (and in N&C page 202).
Here's a hint. To get $|0\rangle$ on the top qubit you'd need either $|00\rangle$ or $|01\rangle$ as the final composite state of the two qubits upon measurement. You can now calculate $\langle 00|U_2U_1|00\rangle$ and $\langle 01|U_2U_1|00\rangle$, where $U_1=H$ and $U_2 = \operatorname{CNOT}$ using the method desribed in the first paragraph of N&C page 202 i.e.
$$\langle y| U_{p(n)}...U_2U_1 |0\rangle = \sum_{x_1,...,x_{p(n)-1}} \langle y |U_{p(n)}|x_{p(n)-1}\rangle \langle x_{p(n)-1}|U_{p_{n-1}}... U_2|x_1\rangle \langle x_1|U_1|0\rangle.$$
Then $|\langle 00|U_2U_1|00\rangle|^2+ |\langle 01|U_2U_1|00\rangle|^2$ gives you the necessary probability. Luckily here we began in the state $|00\rangle$ rather than any arbitrary quantum state $|x\rangle$. Otherwise you'd have to repeat this calculation for all the basis states: $|00\rangle$, $|01\rangle$, $|10\rangle$ and $|11\rangle$.
That is, if $|x\rangle$ was $\alpha|00\rangle+\beta|01\rangle+\gamma|10\rangle + \delta |11\rangle$ rather than $|00\rangle$, the final answer would have been
$|\alpha\langle 00|U_2U_1|00\rangle|^2+ |\alpha\langle 01|U_2U_1|00\rangle|^2 + |\beta\langle 00|U_2U_1|01\rangle|^2+ |\beta\langle 01|U_2U_1|01\rangle|^2 +
|\gamma\langle 00|U_2U_1|10\rangle|^2+ |\gamma\langle 01|U_2U_1|10\rangle|^2 +
|\delta\langle 00|U_2U_1|11\rangle|^2+ |\delta\langle 01|U_2U_1|11\rangle|^2 $
This probability should be $\geq \frac{2}{3}$ for any $x\in\mathcal{A}_{\text{no}}$. Note that you don't have to do any full-fledged matrix multiplication here, as you would have to do in case you wanted to calculate $\operatorname{CNOT}.H|00\rangle$. No unnecessary usage of exponential amount of space!
