Suppose we have some functions $u_0, u_1, u_2: \mathbb{R}_x \times \mathbb{R}_t \rightarrow \mathbb{R}$ that are assumed to be continuous and differentiable in both time and space. I am trying to find the numerical solution for $u_1$ via 4th order Runge Kutta (RK4: https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods) given the equation
$$ \frac{\partial u_1(x,t)}{\partial t} = \sum_{i=0}^2 \sum_{j=0}^2 c_{ij} \thinspace u_i(x,t) \frac{\partial u_j(x,t)}{\partial x} $$
at all grid points $(x_k, t_n)$. Beforehand, I already determined the numerical solutions for $u_0$ and $u_2$ at all these $(x_k, t_n)$ grid points.
The set-up for RK4 is to write
$$\frac{\partial u_1(x,t)}{\partial t} = F(u_1,x,t) = \sum_{i=0}^2 \sum_{j=0}^2 c_{ij} \thinspace u_i(x,t) \frac{\partial u_j(x,t)}{\partial x},$$
where again, the $u_0, u_2$ are "known" (numerically) and $u_1$ is the only "unknown". The $\frac{\partial u_j}{\partial x}$ terms are numeric approximations of the actual derivatives.
Then to determine $\frac{\partial u_1}{\partial t} \big|_{(x,t_{n+1})}$, we have \begin{align*} k_1 = F(u_1(x,t_n),x,t_n) &= \sum_{i=0}^2 \sum_{j=0}^2 c_{ij} \thinspace u_i(x,t_n) \frac{\partial u_j(x,t_n)}{\partial x} \\ &= \sum_{i,j \neq 1} c_{ij}\thinspace u_i(x,t_n) \frac{\partial u_j(x,t_n)}{\partial x} \\ &+ \sum_{j \neq 1} \left(c_{1j} \thinspace u_1(x,t_n) \frac{u_j(x,t_n)}{\partial x} + c_{j1} \thinspace u_j(x,t_n) \frac{u_1(x,t_n)}{\partial x} \right) \\ &+ c_{11} \thinspace u_1(x,t_n) \frac{\partial u_1(x,t_n)}{\partial x}. \end{align*}
What confuses me is the interpretation of $k_2$.
(1) Attempt 1. For known $u_0$ and $u_2$, we input time $t_n + \frac{\Delta_t}{2}$ and only for $u_1$ do we substitute in $u_1(x,t_n) + \tfrac{k_1}{2}.$ That is,
\begin{align*} k_2 &= F(u_1(x,t_n) + \tfrac{k_1}{2},x, t_n + \tfrac{\Delta_t}{2}) \\ &= \sum_{i,j \neq 1} c_{ij}u_i\left(x,t_n + \tfrac{\Delta_t}{2}\right) \frac{\partial u_j\left(x,t_n + \tfrac{\Delta_t}{2}\right)}{\partial x} \\ &+ \sum_{j \neq 1} \left(c_{1j} \thinspace \left( u_1(x,t_n) + \tfrac{k_1}{2} \right) \frac{u_j\left(x,t_n + \tfrac{\Delta_t}{2} \right)}{\partial x} + c_{j1} \thinspace u_j\left(x,t_n + \tfrac{\Delta_t}{2} \right) \frac{\partial}{\partial x} \left( u_1(x,t_n) + \tfrac{k_1}{2} \right) \right) \\ &+ c_{11} \left( u_1(x,t_n) + \tfrac{k_1}{2} \right) \frac{\partial}{\partial x} \left( u_1(x,t_n) + \tfrac{k_1}{2} \right). \end{align*}
The problem with this method is that we don't know values like $ u_0\left(x, t_n + \tfrac{\Delta_t}{2} \right),$ since $u_0$ is only known numerically at the time points $t_n$ and not at half-points $t_n + \tfrac{\Delta_t}{2}$.
(2) Attempt 2. This seems wrong, but we could instead substitute every $u_j(x,t_n)$ with $u_j(x,t_n) + \frac{k_1}{2}$ and evaluate $F$ that way. But then it seems like we are ignoring that $F$ is evaluated at $t_n + \frac{\Delta_t}{2}$, since we are actually still evaluating all the $u_j$ at $t_n$ and then shifting by $k_1.$
EDIT. One strategy of course is to solve for $u_0, u_1, u_2$ simultaneously instead of sequentially. However, the goal here is to use the fact that $u_0$ and $u_2$ were already computed in some black box, so there is no need to re-solve for them. Thank you.
Interpolation by splines seems like your best bet. Since you don't have the values of $u_1$ exactly, you can't get the time derivatives of $u_0$ and $u_2$ exactly, but you can fit a really accurate spline to $u_0$ and $u_2$ and just evalute that. You're going to want to make sure that the spline is more accurate than the RK$4$ scheme so that you aren't introducing error that stays large as the RK$4$ step size goes to $0$. RK$4$ has order $4$, so to be safe you should use a spline of order $5$ or higher, so quartic splines seem like a good starting point.
One problem that can arise here is that at the beginning of the RK$4$ iteration, you may not have enough time points backwards in time to formulate this quartic spline for $u_0$ and $u_2$. One workaround to this might be to solve for $u_0$ and $u_2$ as you have already done, but on a finer time grid for a very short time so you have enough information to approximate the values at the half points to start off the RK$4$ iteration. Once you have done enough time steps of RK$4$, you will have enough points backwards in time to compute these quartic splines and interpolate to get the values at the half points.