@Josolumoh 

I just said that all the moments (including the mean M[1]) are infinite:

M[k] = infinity, so, a "very" closed form!

 Manipulating the integral it is possible to compute this in Maple.
 Unfortunately I encountered some other bugs but they are probably related to the first one. 

Note also that Maple finds a correct (continuous) antiderivative but is not able to compute the limit at infinity.

 Text version:

(1/4)*sqrt(2)*sqrt(1-I)*(1+I)^(3/2)*EllipticE((1/2)*sqrt(2)*sqrt(1+I))+(1/4+1/4*I)*sqrt(2)*sqrt(1-I)*EllipticE((1/2)*sqrt(2)*sqrt(1+I), (1/2)*sqrt(2)-(1/2*I)*sqrt(2))+(-1/2-1/2*I)*EllipticE((1/2)*sqrt(2)-(1/2*I)*sqrt(2), (1/2)*sqrt(2)*sqrt(1+I))-(1/4*I)*sqrt(2)*sqrt(1-I)*EllipticF((1/2)*sqrt(2)*sqrt(1+I), (1/2)*sqrt(2)-(1/2*I)*sqrt(2))+(1/4+1/4*I)*EllipticF((1/2)*sqrt(2)-(1/2*I)*sqrt(2), (1/2)*sqrt(2)*sqrt(1+I))+1/2-(1/8)*sqrt(2)*sqrt(1-I)*(1+I)^(3/2)*EllipticK((1/2)*sqrt(2)*sqrt(1+I));

@tomleslie

It's worth mentioning  that for b[i]=0, the equations are satisfied by any P.

@Markiyan Hirnyk 

All the displayed digits in my comment are correct. It is easy to verify by solving dy/dx=0 with high precision.

@Markiyan Hirnyk 

Actually, with your displayed precision it is

2.0005898444608  at  x = 1.0003097363981.

OK, I have missunderstood the problem. I did not look for the matrix A. 

Actually A is not essential for the content; the author could have started with A without mentioning the method to produce it.

@Markiyan Hirnyk 

Yes. Proof:

@Markiyan Hirnyk  Probably because the author has always LinearAlgebra loaded via an .ini file.
Just load it and it will work.

@Carl Love 

I just obtain

Sol(1);
Error, (in Sol) cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up

plots[odeplot](Sol);
Warning, cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up
Warning, cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up

But the real problem is that using minstep you force dsolve to output something even if the solution does not exist.

@Carl Love 

This approach has a serious flaw. For such an R the ODE has no solution (mathematically). Actually the code will not work unless R(t) is multiplied by a smaller number.
I think that the correct approach is to take R(t) a "continuous random" function e.g.

r:=seq( rand(0.. 0.01)(), i=1..10):

R:=unapply( add(r[i]*cos(10*i*t),i=1..10),t  ):

eq:=diff(a(t),t,t) + a(t) = 3.*(1+R(t))*cos(5.*t):

Sol:= dsolve({eq, a(0)=0, D(a)(0)=1}, a(t), numeric);

plots[odeplot](Sol);

@Carl Love 

The problem is (see the previous answer) that the solution for d is actually d=0, and OP wants d<>0, impossible.

@Josolumoh 

It is definitely not a discrete pdf for an arbitrary (r,b,d).

If the discrete values are {0,1,2,...} then e.g. for [r=3, b=2, d=0.01411053479], F is a pdf.
So, you must choose first the right parameters (and maybe also the discrete values).

@want to be a permanent vegan 

Just add after kitonum's code:

E:=0.5*(y(t)^2+w^2*x(t)^2):
J:=E/w:
plot([E,J], t=0..10, color=[red,blue], thickness=2, legend=['E(t)','I(t)']);

@max125

In 2D input mode you are not forced to use the strict maple syntax. The interface will adjust your input to conform to this syntax before sending it to the maple engine. E.g. a multiplication operator * could be added.

@max125 

You see exactly what you have typed and your input is not parsed.