Question: How to probability density function of a data set of dsolve()?

Dear All,

I would like to plot the probability density function of a state variable obtained from solving differential equations. I have found that there are functions called "PDF" and "KernelDensityPlot" in the Statistics package, but they really confuse me. Could you please point me out? My code is as follows.

Ps. Is it possible to plot the PDF directly from the solution of dsolve() without discretizing the results? 

restart:
with(plots): with(DEtools): with(plottools):with(LinearAlgebra): with(Statistics):

v1:=1: f:=-4: v2:=2.515: omega:=1: epsilon:=0.001: k:=0:
sys:=diff(u1(t),t)=v1*u1(t)-(omega+k*u2(t)^2)*u2(t)-(u1(t)^2+u2(t)^2+3*z(t)^2)*u1(t),
     diff(u2(t),t)=(omega+k*u1(t)^2)*u1(t)+v1*u2(t)-(u1(t)^2+u2(t)^2+3*z(t)^2)*u2(t),
     diff(z(t),t)=z(t)*(-v1+3*u1(t)^2+3*u2(t)^2+z(t)^2)+epsilon*z(t)*(v2+f*z(t)^4):

t_start:=50: t_end:=300: dt:=0.05: fs:=1/dt:

solA:=dsolve({sys, u1(0)=0.6, u2(0)=0.6, z(0)=0.1},
             {u1(t),u2(t),z(t)},
              type=numeric, method=rkf45, maxfun=0,
              output=Array([seq(i,i=t_start..t_end, dt)])):

u1:=solA[2,1][..,2]:
u2:=solA[2,1][..,3]:
z:=solA[2,1][..,4]:

u0:=sqrt~(u1^~2+u2^~2+z^~2):

Phi:=z/~u0:

# How could I plot the probability density of Phi (y-axis) against Phi(x-axis)?

Probability_density_function_plot.mw

Thank you!

Very kind wishes,

Wang Zhe