B is known as a mass matrix. There are many methods which can directly solve with mass matrices. The main ones are implicit methods: implicit Runge-Kutta methods and the multistep BDF methods. Also Rosenbrock methods can handle mass matrices directly, but are not implicit methods.
A quick overview for "why it works" is because at every step these methods have to solve linear equations involving the Jacobian of the RHS. The linear problem is normally of the form $(I - \gamma J)x=b$ where $J$ is the system Jacobian, $\gamma$ and $b$ are some constants, and $x$ is the evolution of the system that must be solved for. Well, it turns out that if you just use the mass matrix $B$ by changing the system to $(B - \gamma J)x=b$, these directly solve the linear problem with the mass matrix, implicitly inverting $B$ in calculations that it would already have been doing (this is because the $I$ comes from the LHS when taking the Jacobian, but with a mass matrix instead that value is just $B$!). Doing this implementation with the non-adaptive 2nd order version of the Rosenbrock method is outlined in Shampine's paper and is not difficult, though you lose a lot of efficiency without adaptive timestepping here.
If you use Julia, DifferentialEquations.jl has access to these methods. Just change
mass_matrix in the construction of the
ODEProblem type. The
Rosenbrock23() method will handle mass matrices like this. In addition, mass matrices can be passed to the implicit RK methods
radau5() which are higher order and tend to be faster if you are looking to solve stiff problems with high accuracy. Technically Sundials'
CVODE_BDF() could solve with a mass matrix, but there's an open issue for setting up the C FFI for doing so (it's high on my priority list. Full disclosure: I am a developer of DifferentialEquations.jl).
If you can handle the performance drop of going to MATLAB instead, you'll find the field is very similar.
ode23s is a 2nd-order Rosenbrock method, is very similar to
Rosenbrock23(), and can use constant mass matrices.
ode15s is very similar to
CVODE_BDF() and can take mass matrices. MATLAB does not have an implicit RK method, though it does have
ode23t which is good for problems where a symplectic integrator is needed and it can handle mass matrices.
I do not see options for passing mass matrices to SciPy's ODE solvers, and have not used the SciPy wrappers. Those wrap some of the same codes as DifferentialEquations.jl so I know that in some cases the mass matrix can be solved by the underlying C/Fortran codes that they use, but it looks like that part of the API is not exposed.