# How can I derive a second order implicit method for this coupled ODE update?

I appreciate this might be an easy question, but I've managed to get myself quite thoroughly confused
I'm trying to solve a system of physics equations that look as follows
$$\frac{\partial \mathbf{E}}{\partial t}=-\alpha \mathbf{j} \\ \frac{\partial \mathbf{j}}{\partial t}=\beta\mathbf{E}+\gamma\mathbf{j}\times\mathbf{B}$$ $$\alpha, \beta, \gamma$$ are scalar constants and $$\mathbf{B}$$ is a constant vector with 3 components. $$\mathbf{E}, \mathbf{j}$$ also both have 3 components.
I would like to solve these equations implicitly. I can easily work out the 1st order implicit Euler method, it ends up with something like $$\begin{pmatrix} \mathbf{E} \\ \mathbf{j} \end{pmatrix}^{n}=\underline{\mathbf{M}}\cdot\begin{pmatrix} \mathbf{E} \\ \mathbf{j} \end{pmatrix}^{n+1}$$ where $$\mathbf{M}$$ is just some 6x6 matrix I can invert through some linear algebra routine. However, I am finding that this implicit Euler method for solving this problem isn't sufficiently accurate and I would like a second order method. I just can't work out how to apply a second order method to this problem.
I've been thinking about using Crank-Nicholson, which for $$\frac{\partial y}{\partial t} = f(t, y)$$ can be written like
$$y_{n+1}=y_n + \frac{h}{2}(k_1 + k_2) \\ k_1 = f(t_n, y_n) \\ k_2 = f(t_n + h, y_n + \frac{h}{2}(k_1 + k_2))$$ But I can't get my head around how you solve the last step for $$k_2$$ for my problem when it isn't so much a function of $$f(t, y)$$ as much as a sort of matrix product like
$$\frac{\partial}{\partial t}\begin{pmatrix} \mathbf{E} \\ \mathbf{j} \\ \end{pmatrix} = \underline{\mathbf{M}}\cdot\begin{pmatrix} \mathbf{E} \\ \mathbf{j} \end{pmatrix}$$ I don't have a good conception of what $$k_1$$ or $$k_2$$ actually look like for my problem.

I think they should look something like $$\mathbf{k_1} = \underline{\mathbf{M}}\cdot\mathbf{y}^n \\ \mathbf{k_2} = \underline{\mathbf{M}}\cdot(\mathbf{y}^n + \frac{h}{2}\underline{\mathbf{M}}\cdot\mathbf{y}^n + \frac{h}{2}\mathbf{k_2})$$ where $$\mathbf{y}$$ is the vector of $$\mathbf{E}, \mathbf{j}$$. But is this remotely correct?

In the linear case you get $$k_1+k_2=2My+\frac h2 M(k_1+k_2),$$ so that isolating the sum you get $$k_1+k_2=2(I-hM/2)^{-1}My.$$ For the step this gives $$y^{n+1}=y^n+\frac h2(k_1+k_2)=y^n+(I-hM/2)^{-1}hMy^n \\ =(I-hM/2)^{-1}[(I-hM/2)+hM]y^n=(I-hM/2)^{-1}(I+hM/2)y^n$$ as one also gets directly from the trapezoidal formula $$\frac{y^{n+1}-y^n}h=\frac{My^n+My^{n+1}}2$$
Practically one could most easily solve for the mean slope $$k_m=\frac{k_1+k_2}2$$, as then $$k_m=\frac{f(t^n,y^n)+f(t^{n+1},y^n+hk_m)}2.$$ If $$J$$ is a sufficiently accurate approximation of the Jacobian around $$(t^n,y^n)$$, then $$f(t^{n+1},y^n+hk_m)-hJk_m$$ will be nearly flat as function of $$k_m$$ and small $$h$$, so that the fixed-point iteration $$k_m^{new}=(I-hJ/2)^{-1}\frac{f(t^n,y^n)+f(t^{n+1},y^n+hk_m^{old})-hJk_m^{old}}2$$ Here "flat" implies a small Lipschitz constant of $$O(h^2)$$ for this iteration. This is usually sufficient to get good enough approximation in 1-3 iterations, as the error in the fixed-point equation decreases by a factor $$O(h^2)$$ in each step. So starting with $$k_m=k_1$$, after one step the error is $$O(h^4)$$ and thus smaller than the $$O(h^3)$$ error for the method step.