bvector

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]);

RK3 step

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