System of ordinary differential equations - time complexity of initial value problem

I am interested in knowing what the time complexity is (in Big-$\mathcal O$ notation) for solving system of $N$ differential equations?

I am using ode15s in MATLAB for this but have no clue about time complexity.

That said, we can easily comment on the time complexity of each step taken by a solver. Supposing we have an IVP $$\dot{x} = f(x,t), \qquad x(0) = x_0$$ with $x \in \mathbb{R}^N$, then for explicit solvers like Matlab's ode45(), the complexity per step is the same as the 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 complexity of $O(N^2)$, since each component of $f(x,t)$ can depend on all components of $x$, but is less in some cases. For example, if $f(x,t)$ is 'sparse' in the sense that each component 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() then a linear system of equations is solved at each step. This involves an $N \times N$ matrix, and the complexity of solving them is $O(N^3)$. This is asymptotically greater than the time required to evaluate $f(x,t)$, so is also the time complexity of the step. Once again though, sparsity can sometimes reduce this for a specific set of ODEs. If 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 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 size to solve a particular IVP, then we could get the time complexity of the whole algorithm by multiplying the above step complexities by $M$, for example $O(MN^3)$ in the implicit case. But usually solvers take steps $\delta t$ of different sizes to achieve a desired accuracy. This makes it impossible to say anything about the overall time complexity without also considering details of the specific IVP. 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 steps while maintaining accuracy. Attempts to analyse this usually consider the 'stability' of the numerical method used on a particular ODE for a particular stepsize, as well as the 'order of accuracy' of the method, and there's a very large literature on the subject.