Trapezoidal method for ODE with a matrix

Hi I know roughly how to use the Trapezoidal method for ODEs but am coming stuck when presents with a system, as follows

$$\bf y'= \begin{bmatrix}\frac{1}{4}x & 1+x^2\\1/2 & -30x\end{bmatrix}+ \begin{bmatrix} cos\pi x \\ sin \pi x \end{bmatrix}$$ also $0<=x<=1$ and $\bf y(0)=\begin{bmatrix} 0\\1 \end{bmatrix}$

So i know what the general case can be expressed , but problems arises when you consider $$\bf W_{i+1} = \frac{2+hl(x_i)}{2-hl(x_{i+1})}W_i+ \frac{g(x_{i+1})+g(x_i)}{2-hl(x_{i+1})}$$ considering $l(x_{i+1})$ is a matrix i am guessing it would have to be a matrix multiplited by another matrix inverse but am struggling to see what that would look like and how to continue on ?

thanks

You just transform $$y_+=y+\frac h2(f(x,y)+f(x_+,y_+))$$ with $f(x,y)=A(x)y+b(x)$ to first get all of $y_+$ on one side and then isolate it $$\left(I-\frac h2 A(x_+)\right)y_+ = \left(I+\frac h2 A(x)\right)y+\frac h2 (b(x)+b(x_+)) \\\iff\\ y_+ = \left(2I - h A(x_+)\right)^{-1}\left[\left(2I + h A(x)\right)y + h (b(x) + b(x_+))\right]$$