The proposed claim $\frac{\partial \exp(DA)}{\partial (DA)} = DA \exp(DA)$ is a category error: it's not even the right kind of object for a matrix derivative. What you need is a Fréchet derivative, a linear operator on matrices, which is not simply a matrix of the same size. See e.g. our matrix calculus course, or this short talk on the nature of the problem.
In particular, the Fréchet derivative $\exp'(X)$ is a linear operator acting on matrices $\delta X$ that gives the first-order change in $\exp(X)$ $$\delta \exp = \exp(X+\delta X) - \exp(X) = \exp'(X)[\delta X] + \mbox{(higher-order)}$$ from a small change $\delta X$ in $X$. Your desired derivative $\frac{\partial \exp(DA) v}{\partial d_k}$ with respect to the $k$-th diagonal entry of $D$ is then $\exp'(DA)[E_k A] v$ where $E_k = \begin{pmatrix} 0 & & & & \\ & 0 & & & \\ & & \ddots & & \\ & & & 1 & \\ & & & & \ddots \end{pmatrix} = \frac{\partial D}{\partial d_k}$ is the diagonal "unit" matrix with a 1 in the $k$-th diagonal.
There are various algorithms to compute a Fréchet derivative operator for a matrix exponential, e.g. Al Mohy et al. (2009) is implemented by JAX in Python and by ChainRules.jl in Julia. As mentioned in the other linked stackoverflow thread, there is also an integral representation that you could implement with a quadrature scheme:
$$
\frac{\partial \exp(DA)v}{\partial d_k} = \frac{\partial \exp(DA)}{\partial d_k} v = \int_0^1 e^{(1-t) DA} E_k A e^{t DA} v \, dt \, .
$$
The basic reason why this doesn't reduce to a much simpler expression is that $E_k A$ doesn't commute with $DA$, so you can't re-order the terms to cancel the $t$.
This formula can be easily checked against a finite-difference approximation, e.g. in Julia for a random $5 \times 5$ case:
julia> using LinearAlgebra, QuadGK
julia> d = randn(5); A = randn(5,5); v = randn(5); D = Diagonal(d); E₁ = Diagonal([1,0,0,0,0]);
julia> f(D,A,v) = exp(DA) v;
julia> (f(D + 1e-8 * E₁, A, v) - f(D, A, v)) / 1e-8 # finite-difference approx.
5-element Vector{Float64}:
-3.6525810709342466
0.04785109530835996
0.012136558424913346
-1.3257712261349752
-0.05547635129055095
julia> quadgk(t -> exp((1-t)DA) * E₁A exp(tDA) * v, 0, 1)[1] # "exact"
5-element Vector{Float64}:
-3.652581057318228
0.047851098313630464
0.012136606329225766
-1.3257711688617824
-0.05547635117441161
(The Al-Mohy algorithm is probably more efficient than quadrature, but the quadrature algorithm is easier to implement yourself to try out. You could also implement the matrix exponential on dual numbers, etcetera. There are famously many ways to compute matrix exponentials, and this is also true for derivatives thereof.)
As mentioned in the other answer, another formulation can be found in Magnus et al (2021) — it actually dates to Mathias (1996). Let me present it in a simpler form. Proposition 2 from that paper shows that the derivative of the matrix exponential $\exp(X)$ with respect to a scalar parameter $x$ is given by the upper-right corner of the exponential of an augmented matrix:
$$
\exp\left(\begin{bmatrix} X & \partial X/\partial x \\ & X\end{bmatrix}\right) = \begin{bmatrix} e^X & \partial (e^X)/\partial x \\ & e^X\end{bmatrix}
$$
which implies that your derivative $\frac{\partial \exp(DA)}{\partial d_k}$ is given by the upper-right corner of
$$
\exp\left(\begin{bmatrix} DA & E_k A \\ & DA\end{bmatrix}\right) = \begin{bmatrix} e^{DA} & \partial (e^{DA})/\partial d_k \\ & e^{DA}\end{bmatrix}
$$
which is easily confirmed in Julia for the $5 \times 5$ example above:
julia> exp([D*A E₁*A; 0I D*A])[1:5,6:10]*v # exact derivative
5-element Vector{Float64}:
-3.652581057318227
0.04785109831363047
0.012136606329225799
-1.3257711688617824
-0.05547635117441162
In contrast, the proposed Jacobian $D A \exp(DA) \mathrm{diag}(Av)$ is incorrect, as is easily verified by comparing to either finite differences or the analytical solution:
julia> (D*A*exp(D*A)*Diagonal(A*v))[:,1] # first col = ∂/∂d₁ ? wrong
5-element Vector{Float64}:
0.41130484171217546
-0.00629410993464498
0.3108523273170691
-2.6517026552656153
0.05596554189606555
