Efficient gate executing the time evolution of a Hamiltonian using Runge-Kutta method

Question

You can find a minimal working example below. In particular, I want to replace the scipy.linalg.expm() matrix exponential by a Runge Kutta time evolution method as this becomes quite slow for a larger system size VQEs.

import numpy as np
import scipy

from qiskit import QuantumRegister
from qiskit.circuit import Gate, QuantumCircuit, ParameterVector



# creates a unitary gate from a list of hamiltonians and parameters
class my_hamiltonian_gate(Gate):
    def __init__(self, num_qubits, hamiltonian_list, theta, label = None):
        self.hamiltonian_list = hamiltonian_list
        super().__init__('my_hamiltonian_gate', num_qubits, theta, label=label)
        
        
    def _define(self):
        qr = QuantumRegister(self.num_qubits)
        qc = QuantumCircuit(qr)
        
        all_qubits = [qr[i] for i in range(self.num_qubits)]
        
        time_evolution_unitary = self.to_matrix()
        qc.unitary(time_evolution_unitary, all_qubits)
        
        self.definition = qc
        
    def to_matrix(self):
        
        hamiltonian = np.zeros(np.shape(self.hamiltonian_list[0]), dtype = np.complex_)
        
        for l in range(len(self.params)):
            curr_param = float(self.params[l])
            curr_hamiltonian = self.hamiltonian_list[l]
            
            hamiltonian += curr_param * curr_hamiltonian
        
        unitary = scipy.linalg.expm(-1j * hamiltonian)

        return unitary



N = 2
qc = QuantumCircuit(N)

param = ParameterVector('a', 1)
param = [param[0]]
hamiltonian = [np.ones((4,4), 'float')]

qc.append(my_hamiltonian_gate(N, hamiltonian, param), list(range(N)))

print(qc)
```