I want to compute the action of matrix exponential on a vector. My matrix $ B$ is of a very large size and it is sparse, e.g., it is written as follows:

size = 350; vector0 = Range[0., 300., 300./(size - 1)]; B0 = NDSolve`FiniteDifferenceDerivative[2, vector0,      "DifferenceOrder" -> 20]["DifferentiationMatrix"];  B = SparseArray@KroneckerProduct[B0, B0]; vector = Flatten@KroneckerProduct[vector0, vector0]; Dimensions[B] 

Based on the documentation and the discussion in https://mathematica.stackexchange.com/questions/83730/possible-method-for-matrixexp/83732?noredirect=1#comment431557_83732, I keep going simply as comes next:

sol2 = MatrixExp[B, vector, Method -> "Krylov"]; // AbsoluteTiming 

This took 31 seconds in my office laptop. It computes the action successfully, BUT it is getting slower and slower when the dimension grows. In the above example, the dimension of $ B$ is $ 122500\times 122500$ , while in my real practice I wish to compute this action for matrices of dimension up to million.

So, I wish to know that is it possible to speed up such a computation by imposing a Compile or by a tolerance or something like that? Since, I am not looking for a result which is correct up to 15 digits! A lower accuracy, maybe up to 6 digits is enough.

Any hints would be welcomed and thanks to Mr. Wizard for encouraging asking this question.