There has definitely been non-zero interest in SDE models over the past few years, with applications spanning physics and biology and epidemiology amongst others, so having a self-contained implementation isn’t a bad idea.
The Euler-Mayorama discretization is very simple but for a stochastic process without much additional structure (like the structure one would get in a Langevin process) it’s the typical tool for which people seem to reach in practice.
One potential issue is with the interface of the function – as most of us don’t use SDEs we’re not particularly well suited to comment on the interface design. In particular should it return a single realization, an ensemble of realizations, moments of the realizations, etc, and then in what format should those outputs be given. I think it would help to have some documentation that more carefully defines the output of the function along with the proposed implementation and a few hypothetical Stan programs using that implementation to help us get a handle on the appropriateness of the interface. Comments from users would also be super helpful.
Finally the biggest concern is the differentiability of the function, which is somewhat related to the interface question above. In general a single realization of a stochastic process won’t be differentiable with respect to the input parameters, hence an ensemble of realizations won’t be either. Conditioning on the stochasticity (i.e. passing in precomputed PRNGs) offers one workaround, as the conditional paths might then be smoothly related to the input parameters. As would returning moments when the discontinuous behavior averages out.
Ultimately because of lack of expertise on these topics within the typical dev community I think we would all benefit from a bit more doc in order to better understand the approach instead of having to deconvolve the method from the code itself. It would certainly help me verify the proposal.