integration of some complicated maximum entropy type expressions

I'm playing with a maximum entropy model for the distribution of momenta 
in the microcanonical ensemble (fixed Energy, N, and V).

My idea is this: based on Jaynes' paper from 1954, I set up the maximum 
entropy distribution over two variables (p the momentum of a specific 
particle, and prest the average momentum of all other particles, I'm 
assuming 1D spatially and an ideal gas with no potential energy just to 
get started)

The constraints are as follows (all expressed nondimensionally):

<(p+prest)/sqrt(2*m*E)> = 0 (on average, no linear momentum)
<(p^2/2m + prest^2/(2*(N-1)*m))/E> = 1 (the total energy is on average E)
<((p^2/2m + prest^2/(2*(N-1)*m))/E - 1)^2> = deltaE (the variance in E 
is a small number)

Eventually I'd like deltaE to go towards 0, making the variance of the 
total energy tend to zero, and then arrive at the joint distribution of 
p and prest.

However, I can't get maxima to calculate the partition function. 
mathematica online will calculate (for example)

integrate(exp(-a*p^2 -b*p^4),p,-inf,inf);

in terms of bessel functions but maxima will not.

This is getting beyond my level of knowledge of analysis techniques. But 
I suspect that if maxima could do this kind of integral it might be able 
to do my problem which is a little more complicated.

---- my exact maxima code follows ----


assume(m > 0);
assume(E > 0);
assume(deltaE > 0);

ppE : exp(-(l[1] + l[2]*(p+prest*(N-1))/sqrt(2*m*E)
     + l[3]*(p^2/2/m)/(E/N) /* this may be an extra constraint I could 
relax?*/
     + l[4]*(p^2/2/m+(prest*(N-1))^2/(2*(N-1)*m) ) /E
     + l[5]*((p^2/2/m + (prest*(N-1))^2/(2*(N-1)*m))-E)^2/E^2
     )
   );

/* this integral is going to fail, I'm sure */
partfcnE:integrate(integrate(ppE,prest,-inf,inf),p,-inf,inf);