How do I determine what $f(y, t)$ is when using Runge-Kutta to solve a PDE? I've been given the formula for RK2 and shown how to use it on an ODE to solve $y'(t) = f(y, t)$, where $f(y,t)$ is explicitly known. 
Equipped only with this knowledge I have been asked to use RK2 on a PDE, where neither side of the equation is known, except that I see that I must substitute  the spacial derivative with a finite difference operator. I am required, for the purpose of this assignment to use a center space finite difference, as shown below.
My work so far: 
Given $\partial_tu(x,t) = -c\partial_xu(x,t)$ is analogous to $y'(t) = f(y,t)$ when we consider all variables besides t as constants:
$$
k_1 = \Delta t f(y_n, t_n) 
\\k_2 = \Delta t f(y_n + k_1, t_n + \Delta t)\\
y_{n+1} = y_n + \frac{1}{2}(k_1 +k_2) + \mathcal{O}(\Delta t^3)
$$
Becomes: 
$$
k_1 = \Delta t f(u_{j+1}, u_{j-1}) = -c\Delta t \frac{u^n_{j+1} - u^n_{j-1}}{2\Delta x} \\k_2 = \Delta t f(u_n +k_1, t_n + \Delta t) = ????? \\
u^{n+1}_j = u^n_j + \frac{1}{2}(k_1 +k_2) + \mathcal{O}(\Delta t^3)
$$
$$
 \text{where }u^n_j \text{ means } u(x_0 + j\Delta x, n\Delta t)
$$
I have no idea how to add $k_1$ to my new "function" of $u_{j+1}$ and $u_{j-1}$ in order to find $k_2$. My textbook (Numerical Recipes [in fortran 77!]) does not review RK2 in combination with PDEs, and neither do my class notes. I have no idea how to functionally put RK2 to work on a PDE without knowing this step. I am not far enough along in the course to make use of linear algebra, though I have seen suggestions to do so. 
This question seems quite similar, but I do not understand the accepted answer:
Runge-Kutta method for PDE
The answer seems to just put finite difference to use. I don't see Runge Kutta used at all. It looks like a finite difference operator is used and that's it.
 A: One approach for performing such time-integration is the method of lines (MOL). Let us introduce the nodal values $u_j(t) = u(x_j,t) = u(x_0 + j\Delta x,t)$. Using centered differencing to approximate the spatial derivative, the PDE gives
\begin{aligned}
\frac{\text d}{\text{d} t} u_j(t) &= -c\, \partial_x u(x_j,t) \\
& \simeq -c\, \frac{u_{j+1}(t) - u_{j-1}(t)}{2\, \Delta x} \, .
\end{aligned}
The system of ordinary differential equations
\begin{aligned}
\frac{\text d}{\text{d} t} {\bf u}(t) &= -\frac{c}{2\, \Delta x}
\left(
\begin{array}{ccccc}
 0 & 1 & &\\
-1 & \ddots & \ddots &\\
 & \ddots & \ddots & 1\\
 & & -1 & 0
\end{array}
\right)
\, {\bf u}(t) - \frac{c}{2\, \Delta x}
\left(
\begin{array}{c}
-u_0(t)\\
0\\
\vdots\\
0\\
u_N(t)
\end{array}
\right) ,\\
&= {\bf f}({\bf u}(t)) \, ,
\end{aligned}
made by the vector ${\bf u} = (u_1,\dots ,u_{N-1})^\top$ of nodal values can then be integrated in time explicitly using Runge-Kutta methods.
In particular, if the forward Euler (RK1) method ${\bf u}^{n+1} = {\bf u}^{n} + \Delta t\, {\bf f}({\bf u}^n)$ is used for time-integration, then an unstable scheme
$$
u_j^{n+1} = u_j^n -c\, \frac{\Delta t}{2\,\Delta x} (u_{j+1}^n - u_{j-1}^n)\, , \qquad 1 \leq j \leq N-1
$$
is obtained. A slight modification of this method gives the stable Lax-Friedrichs scheme.
If the proposed improved Euler (RK2) method is used for time integration, one has
\begin{aligned}
\tilde u_j^{n+1} &= u_j^n + \Delta t\, k_1 \\
&= u_j^n -c\, \frac{\Delta t}{2\, \Delta x} (u_{j+1}^n - u_{j-1}^n) \, ,\\
u_j^{n+1} &= u_j^n + \frac{1}{2} (k_1+k_2)\\
&= u_j^n -c\, \frac{\Delta t}{4\, \Delta x} (u_{j+1}^n - u_{j-1}^n + \tilde u_{j+1}^{n+1} - \tilde u_{j-1}^{n+1}) \, .
\end{aligned}
Finally, for $2 \leq j \leq N-2$,
$$
u_j^{n+1} = u_j^n -c\, \frac{\Delta t}{2\, \Delta x} (u_{j+1}^n - u_{j-1}^n) + c^2\, \frac{\Delta t^2}{8\, \Delta x^2}(u_{j+2}^n - 2u_{j}^n + u_{j-2}^n) \, ,
$$
which looks quite similar to the Lax-Wendroff scheme, but may have different stability properties.
