I'm struggling to try and put my idea of what I have for this problem into Python, I'm stuck on trying to put the bvector(x) function to give me my required output.
This encodes the inhomogeneity in the linear system of differential equations $$ y'(x)=A\,y(x)+b(x) $$ Thus $b(x)$ is a vector of the same dimension as $y$. As the vector addition in matlab and python has a mode of adding a scalar, adding it to all components, one can simplify the current case of a zero vector to just returning the scalar $0$.
bvector = lambda x: 0
If this vector is non-trivial, you would have to return a proper vector, for instance using
numpy.array as vector type. In dimension two this could look like
bvector = lambda x: np.array([x+1,x**2]);
If the current state is
xn,yn then the step of size
h is implemented as just repeating the (corrected) formulas,
rk3step(xn,yn,h): y1 = yn+h*(A.dot(yn)+bvector(xn)); y2 = 0.75*yn+0.25*y1+0.25*h*(A.dot(y1)+bvector(xn+h)); return (yn+2*y2+2*h*(A.dot(y2)+bvector(xn+0.5*h)))/3;
A possible loop
This you now can embed into a list construction loop
def rk3(A, bvector, y0, interval, N): def rk3step: ... x = np.linspace(*interval,N+1); y = [y0] for n in range(N): y.append(rk3step(x[n],y[n],x[n+1]-x[n])); return x,np.asarray(y).T