While experimenting with some linear algebra and some integrals, I believe I may have found one way of evaluating the general antiderivatives of $e^{ax}\sinh\left(bx\right)$ and $e^{ax}\cosh\left(bx\right)$. Hopefully, this is not a duplicate as I have used Approach0 to try to find similar methods but could not find any. I apologize in advance if this is a lot to read, so I'll try my best to be as concise as possible.
(Proposal) Evaluate $\displaystyle\int e^{ax}\cosh\left(bx\right)dx$ and $\displaystyle\int e^{ax}\sinh\left(bx\right)dx$ where $a,b \in \mathbb{R}$ and $a,b$ are both not $0$.
(Attempt) Consider
$$\operatorname{span}\left\{e^{ax}\cosh\left(bx\right), e^{ax}\sinh\left(bx\right)\right\} := \operatorname{span}\left\{v_1,v_2\right\}$$
and the differentiation operator
$$D(f) := f'.$$
To prove the span of the two vectors is invariant under $D$, observe that
$$D(v_1) = \frac{d}{dx}e^{ax}\cosh\left(bx\right) = ae^{ax}\cosh\left(bx\right)+be^{ax}\sinh\left(bx\right)$$
and
$$D(v_2) = \frac{d}{dx}e^{ax}\sinh\left(bx\right) = be^{ax}\cosh\left(bx\right)+ae^{ax}\sinh\left(bx\right).$$
Thus, applying $D$ to the two basis vectors still yields a linear combination of those basis vectors, so we can find a matrix for $D$ by gathering the coefficients. Let
$$A := \begin{bmatrix} a & b \\ b & a \\ \end{bmatrix}.$$
Then
$$D\begin{bmatrix} v_1 \\ v_2 \\ \end{bmatrix} = A\begin{bmatrix} v_1 \\ v_2 \\ \end{bmatrix}.$$
To antidifferentiate, we find the inverse of $D$, call it $D^{-1}$, so $D^{-1} = A^{-1}.$ Suppose $\overrightarrow{C} = \begin{bmatrix} c_1 \\ c_2 \\ \end{bmatrix} \in \mathbb{R^2}.$ Then
$$\displaystyle \int \begin{bmatrix} v_1 \\ v_2 \\ \end{bmatrix}dx = A^{-1}\begin{bmatrix} v_1 \\ v_2 \\ \end{bmatrix} + \overrightarrow{C} = \dfrac{1}{\operatorname{det}(A)}\begin{bmatrix} a & -b \\ -b & a \\ \end{bmatrix}\begin{bmatrix} v_1 \\ v_2 \\ \end{bmatrix} + \overrightarrow{C}.$$
Simplifying gives us
$$\displaystyle \int \begin{bmatrix} v_1 \\ v_2 \\ \end{bmatrix}dx = \frac{1}{a^2-b^2}\begin{bmatrix} av_1-bv_2+c_1 \\ -bv_1+av_2+c_2 \\ \end{bmatrix}.$$
Multiplying $\dfrac{1}{a^2-b^2}$ into the vector provides us the general antiderivative of $v_1$ on the top and the general antiderivative of $v_2$ on the bottom.
(Question) Is this correct? If so, how can we apply this method to other integrals out there?
Any answer/comment shedding some light is appreciated.
(Thoughts) It seems like this method only works if we can solve $D(\overrightarrow{v}) = A(\overrightarrow{v})$, where $D, A$ are the same as in my solution and $\overrightarrow{v}$ is an element of the span of two elements. This leads me to believe I was somewhat lucky in getting the right antiderivatives. For context, my background in algebra consists of an undergraduate, upper-division linear algebra course, and an abstract algebra course. If I missed anything trivial, please let me know.