I'm going to assume that you're specifically talking about Initial value Problems (IVPs) for Ordinary Differential Equations (ODEs), since that's what ode15s() is used for.
This is quite a complicated question to answer, because good ODE solvers have stepsize control procedures that automatically take smaller or larger steps to achieve a desired accuracy. As a result the actual amount of time taken to solve the system of ODEs can vary a lot depending on the specific system (and sometimes, initial conditions) considered. Also, the big-O notation hides constants which are often important when comparing different methods.
That said, we can easily comment on the complexity of each step of an ODE solver. Supposing we have an ODE $$ \dot{x} = f(x,t) $$ with $ x \in \mathbb{R}^N $, then if you use an explicit solver like Matlab's ode45(), the time complexity of each step is the same as the time complexity of evaluating $f(x,t)$. This is because explicit solvers evaluate $f(x,t)$ a fixed number of times per step. Typically this means a time complexity per step of $O(N^2)$, since each component of $f(x,t)$ can depend on all components of $x$, but is possibly less in some cases. For example, if $f(x,t)$ is 'sparse' in the sense that each component of the result only depends on at most $m$ components of $x$, then the time complexity per step drops to $O(mN)$.
On the other hand, if you use an implicit solver like ode15s(), as in your case, a linear system of equations is solved at each timestep. This system of equations involves an $N \times N$ matrix, so the complexity of solving them is $O(N^3)$. This is asymptotically greater than the time required to evaluate $f(x,t)$, and so is also the time complexity of the whole step. Once again though, sparsity of the linear system can sometimes reduce this for a specific set of ODEs. Again considering the case that each component of $f(x,t)$ only depends on at most $m$ components of $x$, we can often construct the linear system such that it involves a banded matrix of bandwidth $m$, and then the time complexity per step reduces to $O(m^2N)$.
Now here's where it gets complicated. If we were to fix a stepsize $\delta t$ and take $M$ steps of that fixed stepsize to solve a particular IVP, then we would get the time complexity of the whole solve procedure by multiplying the above step complexities by $M$, for example $O(MN^3)$ in the implicit case. But usually solvers take steps of different sizes $\delta t$ to achieve a desired accuracy. This makes it impossible to say anything about the overall time complexity without also considering some details of the specific problem being solved. And despite the fact that implicit methods like ode15s() have a worse time complexity per step than explicit methods like ode45(), for some problems (called stiff problems) we find that implicit methods outperform explicit methods as they can take larger stepsizes. Attempts to analyse this usually consider the 'stability' of the numerical method used on a particular problem for a particular stepsize, which can be evaluated by linearising $f(x,t)$, as well as the 'order of accuracy' of the method, and there's a very large literature on the subject.
In practice, the usual procedure is to first try an explicit solver like ode45() (which is the Dormand-Prince method) and then if the problem is taking a long time to solve due to stiffness, try an implicit method like ode15s(). Since there's a huge number of different methods and software for solving IVPs, the only way to find out the best method for a particular problem is benchmarking.