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$.
Best Answer
Have you tried this online MMA Web Based tool to do a 3D Plot of Vector Fields?
You might also try have a look here and see if any of the free SW listed works for you (GNU Plot, JavaView, Visual Data...).
Have fun!
Regards