It's not terribly pretty, but here's a start in sage:
p = [] G = Graphics() xw = 4 # maximum/window "radius" for x yw = 10 # maximum/window "radius" for y var('x,m') assume(x > 0) assume(m > 0) y = e^(m*x^2/2) p.append(y*integrate(1/y,(x,0,x))) G += plot(p[-1].subs(m=1),-xw,xw) p.append((y-1)/m) G += plot(p[-1].subs(m=1),-xw,xw,color='green') for k in range(3): p.append(p[k]*(k+1)/m-x^(k+1)/m) G += plot(p[-1].subs(m=1),-xw,xw,color=(1,k/3,k/3)) G.show(ymin=-yw,ymax=yw,figsize=4)

Here's two versions of the above code, made into a function to play around with. This first version requires the $x$ and $y$ limits. It also shows you the functions it calculated.
def myplot(m,n,x_max,y_max): p = [] G = Graphics() var('x') assume(x > 0) y = e^(m*x^2/2) p.append(y*integrate(1/y,(x,0,x))) G += plot(p[-1],-x_max,x_max) p.append((y-1)/m) G += plot(p[-1],-x_max,x_max,color='green') for k in range(n): p.append(p[k]*(k+1)/m-x^(k+1)/m) G += plot(p[-1],-x_max,x_max,color=(1,k/n,k/n)) G.show(ymin=-y_max,ymax=y_max) return p
This second version does not require $y$ limits, but scales automatically in the $y$ dimension.
def my_plot(m,n,x_max): p = [] G = Graphics() var('x') assume(x > 0) y = e^(m*x^2/2) p.append(y*integrate(1/y,(x,0,x))) G += plot(p[-1],-x_max,x_max) p.append((y-1)/m) G += plot(p[-1],-x_max,x_max,color='green') for k in range(n): z = p[k]*(k+1)/m-x^(k+1)/m p.append(z) G += plot(z,-x_max,x_max,color=(1,k/n,k/n)) G.show()
If you can, try playing around with these and let me know what you find. There are a few potential problems. (1) is it calculating a formula that works for $x<0$? (2) I draw $P_0$ in blue, $P_1$ in green, and the iterates in successively lighter shades of red. However, I didn't seem to see as many shades of red as I expected. It should draw $n$ iterates, $n$ being the second function argument. (3) There are some runtime errors but it does produce a plot. For example:
myplot(-1,5,6,6)

Hope this helps for now! I'll try to check back within 24 hours.
Upon reflection/recollection, I noticed and was going to ask whether you are aware that $P_{n}(x)$ is the particular solution to the differential equation $y'-mxy=x^n$ with initial condition $(0,0)$, but I see now that you are!
You can also characterize the qualitative behavior for a fixed $m$ from the vector field plot, which easily justifies your intuition that $y$ is bounded when $m<0$ and grows rapidly when $m>0$.