First of all, don't even consider "optimizing" before you're using the right integration method. Chucking more computers at the problem may sound like the easiest way to solve it, but in reality it will be more difficult than you'd think (due to implicit BLAS multithreading explained later). If you're using explicit Euler, you can do a whole lot better just by changing the method. Proper adaptive time stepping (possibly implicit depending on the stiffness) algorithms can be orders of magnitude more efficient than explicit Euler. Don't even consider parallelizing and keeping explicit Euler unless you have a very good reason (i.e. for hyperbolic PDEs so you're concerned about SSP, but in that case just use an SSP-optimized integrator... you're not going to have a good reason).
Next, your derivative calculation function is by far the most expensive part of the calculation. First of all, you should make sure that it's optimized. Language won't matter much here since all of the time will be in the BLAS kernel for matrix multiplication, but make sure it's non-allocating via in-place operations, using BLAS or SIMD'd loops as much as possible, etc.
Once that's all good, consider whether you still need to parallelize. With the right ODE methods and with a decently performant derivative implementation, 10^5 can be reasonably speedy. If you ended up using an implicit method, you'll want to do things like swap out the linear solver to some iterative solver, and that will give you another nice speedup.
But since this is a post about parallelism, I'll keep going, assuming you want to do something extreme.
How to parallelize this: the details to know
Since $A$ and $D$ are matrices, if you're using a language with a standard BLAS setup then you're already parallelizing via multithreading. Just check your CPU usage (ex:
htop on Linux) and you'll see that large enough matrix multiplications with multithread in Python, R, Julia, MATLAB, etc. So, getting more performance out isn't as simple as just parallelizing it since if you're on a single computer it's likely you're using all of your cores optimally already.
Parallel in time algorithms like parareal only make sense when the ODE is cheap since there's no other available parallelism. Parallel in time algorithms are not efficient, but they try to make up for it by allowing parallelism. The overhead cutoff can be quite high for these algorithms, I've read that 32 cores is a good low estimate. These are definitely not suitable to the way you've described your problem.
However, there's two other levels of parallelism. If your solving is expensive, then solving at multiple systems at the same time in parallel will be almost optimal speedup. i.e. do a parallel for loop or map that solves 16 systems at a time in a 16 core machine and it will be about 16x speedup since the data transfer times are almost 0 compared to the full integration. For example, if you're using Julia, you can write this up yourself quite easily, or use what's built into DifferentialEquations.jl, or find the right resource for whatever other language. Remember, BLAS is already multithreading (and most likely those matrix multiplications are the most expensive operation in the loop) so this only makes sense if you can farm out to other computers/nodes. Julia's native parallelism does that, Python's multiprocessing can do it, you can use MPI in C/Fortran, etc. If you have a big cluster and want lots of gains for little work, this is a way to solve your problem.
Finally, if you have an optimized derivative function with most of the time in matrix operations but don't have extra nodes hanging around (or if you do, you can still apply this), another level of parallelism you can exploit is within-method. You're already doing within-level parallelism via the multithreaded BLAS operations, but you can do better. GPUs are essentially parallel linear algebra machines, so if you have one hanging around then you can just perform the matrix multiplications on then it's a good idea. Setting it up isn't hard since things like CuBLAS are just BLAS on GPUs, so you just have to call it. However, I think you'll have memory issues if you do it naively, so this could be difficult.
Stop using explicit Euler: the first thing you should do is really optimize your code and your method. If that's not enough, solve the system on different nodes at the same time. If that's not enough, accelerate the linear algebra with GPUs/Xeon Phis/etc. since the limiting step is the matrix operations.
For reference, here's a notebook explaining the optimizing step and here's a blog post explaining within-method GPU parallelism.