Need help with a Crank Nicholson Method example problem. I have an exam coming up and the professor released the sample test containing a Crank Nicolson question.  I was out of town for those two lectures, so I missed the information.  Even though I have acquired the notes, the professor didn't do an example problem, which is the best way I learn a new method.  The problem is as follows:

Apply the Crank-Nicolson method to the problem
  \begin{align}
\frac{∂u}{∂t}&=\frac{∂^2u}{∂x^2},& &0\le x\le 1,\quad t>0\\
u&=x(x-1),& &0<x<1,\quad t=0\\
\frac{∂u}{∂x}-u&=e^{-t},& &x=0,\quad t>0\\
\frac{∂u}{∂x}&=3& &x=1,\quad t>0
\end{align}
  Let $h=\frac12$, $k=\frac12$ and write the implicit formula for the $(n+1)$th row. Invert the matrix and find an explicit form.

I've been trying to catch up on how to do this method (I watched this video, and I've read a few written example problems), but I'm still not really grasping the process.  I know I need to create the grid and from that derive matricies, but I don't really understand how.  I really need a solved example in order to understand algorithms such as these.
I was hoping someone could help me solve the above sample problem so that I can understand the Crank Nicolson method for my upcoming exam.
Basically the only thing I've done right (I think...) is set up a 2x2 grid that goes from 0 to 1 in 1/2 increments.
 A: Crank-Nicolson for the PDE $\frac{∂u}{∂t}=L(t,u)$ plus boundary conditions first performs a discretization in the space dimension(s) forming a vector of sample point values $\vec u(t)=\{u(t,x_k)\}_{k=0,…,N}$, $x_k=x_0+k·Δx$, and transforming the operator $L$ plus the boundary conditions into a function $f(t,\vec u)$ of this vector approximating them.
Then it uses the trapezoidal implicit Runge-Kutta or multi-step method
$$
\vec u(t+Δt)= \vec u(t)+\frac{Δt}2·\Bigl(f(t,\vec u(t))+f(t+Δt,\vec u(t+Δt))\Bigr)+O((Δt)^3)
$$
or
$$
2\vec u(t+Δt)-Δt·f(t+Δt,\vec u(t+Δt))=2\vec u(t)+Δt·f(t,\vec u(t))
$$
to compute $\vec u(t+Δt)$ from $\vec u(t)$.

The building blocks for the construction of $f$ from $L$ are
he second order central difference quotient
$$
\frac{∂^2u}{∂x^2}(t,x)=\frac{u(t,x-Δx)-2u(t,x)+u(t,x+Δx)}{(Δx)^2}+O((Δx)^2)
$$
and the first order central difference
$$
\frac{∂u}{∂x}(t,x)=\frac{u(t,x+Δx)-u(t,x-Δx)}{2·Δx}+O((Δx)^2)
$$
to approximate (virtual) values for $u(t,x_{-1})$ and $u(t,x_{N+1})$, 
or back- and forward differentiation formulas to incorporate the boundary values and boundary conditions. This combines, as one possibility, to the function
$$
f(t,\vec u)=(f_k(t,\vec u))_{k=0,...,N},
$$
with components
$$
f_k(t,y)=\begin{cases}
\dfrac{u_{k-1}-2u_k+u_{k+1}}{(Δx)^2}&k=1,…,N-1\\[1em]
\dfrac{\bigl(u_1-2Δx·\bigl(u_0+e^{-t}\bigr)\bigr)-2u_0+u_{1}}{(Δx)^2}&k=0\\[1em]
\dfrac{u_{N-1}-2u_N+\bigl(u_{N-1}+2Δx·3\bigr)}{(Δx)^2}&k=N\\
\end{cases}
$$

Ultimately, the operator $L$ and the function $f$ are linear in $u$ resp. $\vec u$, at least in this example. Thus $f$ can be written in matrix form $$f(t,\vec u)=A(t)·\vec u+b(t).$$
For the time evolution this means 
$$
$$
\Bigl(I-\tfrac12Δt·A(t+Δt)\Bigr)·\vec u(t+Δt)=\Bigl(I+\tfrac12Δt·A(t)\Bigr)·\vec u(t)+\tfrac12Δt·(b(t)+b(t+Δt)).
$$
The matrices are tri-band, so it makes in general no sense to write down an explicit inverse. However, in the toy example with 3 sample points this is still, if barely, feasible and seems to be required in the task.
