Recursive solutions to linear ODE.
When finding the solutions to the simple ODE
$$ y'- mxy= x^n \text{ ; } y(0) = 0$$
I found the following:
Let $P_n$ be the particular solution for each integer exponent $n$. Then if we define
$$P_0(x) = \exp\left(\displaystyle \frac{mx^2}{2}\right)\int_0^x \exp\left(-\displaystyle \frac{mt^2}{2}\right)dt$$
$$P_1(x) =\frac{1}{m}\left\{ \exp\left(\displaystyle \frac{mx^2}{2}\right)-1 \right\}$$
All other solutions are given by the following recursion
$$P(x)_{n+1} = \frac{n}{m}\left\{ P_{n-1}(x)-\frac{x^{n}}{n} \right\} $$
$$P'(x)_{n+1} = nxP(x)_{n-1}$$
Is there any theory on such functions? (I'm mostly interested in the last one, which is very similar to Bernoulli's:
$$B'_{n+1}(x) = (n+1)B_{n}(x)$$
ADD: The particular solutions are
$${P_{2n + 1}} = \frac{{\left( {2n} \right)!!}}{{{m^n}}}{P_1} - \sum\limits_{k = 1}^n {\frac{{\left( {2n} \right)!!}}{{\left( {2k} \right)!!}}} \frac{{{x^{2k}}}}{{{m^{n - k + 1}}}}$$
$${P_{2n}} = \frac{{\left( {2n - 1} \right)!!}}{{{m^n}}}{P_0} - \sum\limits_{k = 1}^n {\frac{{\left( {2n - 1} \right)!!}}{{\left( {2k - 1} \right)!!}}} \frac{{{x^{2k - 1}}}}{{{m^{n - k + 1}}}}$$
So maybe the importance will be in the polynomials (the sums).
Solution in terms of the incomplete gamma function
The behavior of the solutions to this ODE are somewhat obscured by its representation. Let $t = m x^2/2$ and $u_a = 2(m/2)^a P_{2a-1}$, where $a = (n+1)/2$. The ODE takes the form $$u_a'(t) - u_a(t) = t^{a-1},$$ with the boundary condition $u_a(0) = 0$. The Green's function for the operator $d/dt - 1$ satisfying the boundary condition is $G(t,t') = \Theta(t-t') e^{t-t'}$. Therefore, $$u_a(t) = e^t \int_0^t d z\, e^{-z} z^{a-1}.$$ The integral is the lower incomplete gamma function, so \begin{equation*} u_a(t) = e^t \gamma(a,t). \tag{1} \end{equation*} In terms of $P$ and $x$ we find $$\begin{eqnarray*} P_n(x) &=& \textstyle\frac{1}{2} \left(\frac{2}{m}\right)^{\frac{n+1}{2}} u_{\frac{n+1}{2}}\left(\frac{m x^2}{2}\right) \\ &=& \textstyle\frac{1}{2} \left(\frac{2}{m}\right)^{\frac{n+1}{2}} \exp\left(\frac{m x^2}{2}\right) \gamma\left(\frac{n+1}{2}, \frac{m x^2}{2}\right). \end{eqnarray*}$$
The interesting properties of $P_n$ are inherited from $\gamma$. For example, $$\begin{eqnarray*} u_a(t) &=& a^{-1} t^a M(1,a+1,t) \\ &=& \Gamma(a) t^a \sum_{k=0}^\infty \frac{t^k}{\Gamma(a+k+1)} \\ &=& \frac{t^a}{\displaystyle a-\frac{a t}{\displaystyle a+1+\frac{t}{\displaystyle a+2-\ldots}}}. \end{eqnarray*}$$ Above $M(a,b,z) = {}_1 F_1(a;b;z)$ is Kummer's confluent hypergeometric function. Have a look at Wikipedia and the DLMF on the incomplete gamma function and related functions.
In the question statement only integer and half-integer values of $a$ were discussed. Equation (1) above defines the function for all $\mathrm{Re}\, a >0$. In fact we can analytically continue to most of the complex plane, in which case $u_a(t)$ is meromorphic with simple poles at $a = 0, -1, -2, \ldots$. For convenience, $$\begin{eqnarray*} u_{1/2}(t) &=& \sqrt{\pi} \mathrm{erf}(\sqrt{t}) e^t \\ u_1(t) &=& -1 + e^t. \end{eqnarray*}$$ These are $P_0$ and $P_1$ in terms of $u$ and $t$.
Recurrence relation and derivative
The recurrence relation given in the statement of the question can be thought of as a consequence of that for $\gamma$, $$\gamma(a,t) = (a-1)\gamma(a-1,t) - t^{a-1}e^{-t}.$$ (This can be found by integration by parts.) Thus, we find $$u_{a+1} = a u_a - t^a.$$ The identity of interest takes the form $$\frac{d u_{a+1}}{d t} = a u_a.$$ Notice how the original representation obscures this relation. (Note that the formula for the Bernoulli polynomials should read $B_{n+1}' = (n+1)B_n$.) The boundary conditions for these two functions are completely different. The Bernoulli polynomials satisfy $B_n(0) = B_n$ whereas we have $u_a(0) = 0$. The sequence of functions $u_a$ is not an Appell sequence. We can't index the $u$ in any way such that $u_0$ is constant. Of course, the $u$ are not polynomials.
A related Appell sequence
However, if we change the boundary conditions and consider $v_a(t) = u_{a+1}(t) + c_a e^t$ (so that $v_a(0) = c_a$) and if $c_a = -a!$ the sequence is Appell, $$v'_{a+1} = (a+1)v_a, \hspace{5ex} v_0(t) = -1.$$ (Note that $v_{a-1}$ satisfies the differential equation $v' - v = t^{a-1}$.) These polynomials are not very interesting, $$\begin{eqnarray*} v_a &=& -e^t \Gamma(a+1,t) \\ &=& -a! e_a(t), \end{eqnarray*}$$ where $e_a(t) = \sum_{k=0}^a t^k/k!$ is the partial exponential sum. (Here we have assumed $a$ is integral. For half-integral $a$ this does not work.)
Polynomials associated with original boundary conditions
We now consider the polynomials associated with the original boundary condition. For integral $a$ we let $q_n = u_{n+1} - n!u_1$. Then $$q_n = -n! \sum_{k=1}^n \frac{t^k}{k!}.$$ These are just the polynomials $v$ found above without the constant term. Note that $\frac{1}{2} (\frac{2}{m})^{n+1} q_n(\frac{m x^2}{2}) = - \sum_{k=1}^n \frac{(2n)!!}{(2k)!!} \frac{x^{2k}}{m^{n-k+1}}$, since $n! = (2n)!!/2^n$.
For half-integral $a$ we let $r_n = u_{n+1/2} - \frac{\Gamma(n+1/2)}{\Gamma(1/2)} u_{1/2}$. Then $$r_n = -\Gamma(n+1/2) \sum_{k=1}^n \frac{t^{k-1/2}}{\Gamma(k+1/2)}.$$ Note that $\frac{1}{2} (\frac{2}{m})^{n+1/2} r_n(\frac{m x^2}{2}) = - \sum_{k=1}^n \frac{(2n-1)!!}{(2k-1)!!} \frac{x^{2k-1}}{m^{n-k+1}}$, since $\Gamma(n+1/2) = (2n-1)!!\sqrt{\pi}/2^n$.
The sums for $q$ and $r$, though finite, are fairly ordinary. Perhaps surprisingly so. Their coefficients aren't related to Bernoulli numbers or any other interesting sequences. It doesn't appear useful to study these further. In fact, it is somewhat unnatural to separate the solution into two parts that fail to solve the same differential equation.
Further exploration
To explore further the functions $P_n$ it is probably more fruitful to study the properties of the incomplete gamma function and related special functions in more detail. The DLMF will be very useful in this regard. Here are two examples of the sorts of relations one can find. (It is probably best to stick with $u$ and $t$, but one can always go to $P$ and $x$ using the dictionary given at the top of this document.)
Using the properties of $\gamma$ we find $$u_{a+n}(t) = (a)_n u_a(t) - t^a \sum_{k=0}^{n-1} \frac{\Gamma(a+n)}{\Gamma(a+k+1)} t^k,$$ where $(a)_n = \Gamma(a+n)/\Gamma(a)$ is Pochhammer's symbol. The recurrence relation is a special case of this formula when $n=1$. Lastly, a somewhat amazing formula, $$u_a(\lambda t) = \lambda^a \sum_{k=0}^\infty u_{a+k}(t) \frac{(1-\lambda)^k}{k!}.$$ There is much more to learn about your functions by exploiting the properties of $\gamma$. Cheers!