options = odeset('mass',BigM,'MassSingular','no','BDF','on','MaxOrder',2,'Jacobian',...
@(t,u) MyJacobian(t,u,MM_params, RHS_params, mask));
[time,u_sol] = ode15s(@(t,u) TGetRHS(t,u, MM_params, RHS_params, mask), 0:0.4:tspan, U_init, options);
I have a project in which the main bulk of the work is solving a huge ODE coupled system of size 
of the form: where 
The RHS is partitioned into the sum of a "fast scale" and "slow" scale component. This makes the system very stiff, hence I am using ode15s in MATLAB and I tried to make things as sparse as possible.
However, my code takes several hours to run, and I have been thinking about how to use parallel or gpu computing to speed things up. I have checked MATLAB documentation to see if Ode15s can be solved on GPU, but I havent found any good answer.
I will appreciate any suggestion on how to make my code run faster in MATLAB or using other external software.
