I don't know why this is unexpected.  Your polynomial is a quartic in z = Z^2 with general coefficients except for no linear term (i.e. you could get any quartic  c4*z^4 +c3*z^3+c2*z^2+c0 by appropriate choice of A,B,C,D,L).  A suitable translation z -> z - k will give you a completely general quartic.  The formula for the solution of the general quartic in radicals is extremely complicated.

The solution of the N x N matrix Riccati system

P' = A+B.P+P.B^%T+P.C.P

with initial value P(0) = P0 is P(t) = U(t) . V(t)^(-1) where

[ U ]         [ B    A    ]  [ P0 ]
[ V ] = exp(t [ -C  -B^%T ]) [ Id ]

(Id being the N x N identity matrix, and ^%T meaning transpose).  Unfortunately if you try
to compute this in closed form it will probably be extremely complicated.  Even using floating
point, although somewhat better, the formulas will be quite complicated.  Here is an example.

>  with(LinearAlgebra): N:= 4;
   A:= RandomMatrix(N,N,generator=rand(-2..2),outputoptions=[shape=symmetric]);
   B:= RandomMatrix(N,N,generator=rand(-2..2));
   C:= RandomMatrix(N,N,generator=rand(-2..2),outputoptions=[shape=symmetric]);
   P0:= RandomMatrix(N,N,generator=rand(-10..10));
   H:= <<B | A>,<-C | -B^%T>>;
   Digits:= 14;
   M:= MatrixExponential(evalf(H),t);
    

 You may notice that although the coefficients should be all real, there are some small imaginary parts here due to roundoff error.
To fix that:

> M:= map(evalc @ Re, M);
  USVS:= M . <P0 , IdentityMatrix(N)>:
  US:=USVS[1..N,1..N]:
  VS:=USVS[N+1..2*N,1..N]:
  PS:= map(combine, US . Adjoint(VS))/combine(Determinant(VS)):
  length(PS);
    

In this case the length I get is 847012: rather big!

This particular solution, it turns out, has a singularity at about t = 0.0232.

> eval(PS, t = 0.02);

        [-5.0522133792393 , 9.4687046624662 , 60.613559013266 , -13.976165888432]

        [-27.836390227848 , 16.654194978316 , 30.463674741680 , -13.433874731242]

        [-7.4364059647309 , -1.3214480205026 , 2.2174291183686 ,  -3.2024126939127]

        [-49.268210843687 , 23.127583272672 , 83.227607337801 , -17.911747659531]

For comparison, here's a purely numerical solution.

> P:= Matrix(4,4,(i,j) -> p[i,j](t));
  sys:= convert(map(diff,P,t) - (A + B . P + P . B^%T + P . C . P), set);
  ics:= {seq(seq(p[i,j](0)=P0[i,j],i=1..N),j=1..N)};
  PSN:= dsolve(sys union ics, numeric);
  eval(P, PSN(0.02));


        [-5.05221376779952180 ,  9.46870493177853412 , 60.6135598639016494 , -13.9761661205448836]

        [-27.8363904078385610 ,  16.6541950981618890 , 30.4636751360981428 , -13.4338748339841381]

        [-7.43640599979502337 , -1.32144799151080372 , 2.21742917653212634 , -3.20241269643678494]

        [-49.2682114316373117 ,  23.1275836655751803 , 83.2276086063742753 , -17.9117479934783966]

The int function tries to do an exact symbolic integration.  Since your expression is extremely complicated, it's unlikely to succeed.  To get a numerical answer, use Int rather than int and then use evalf on the result.  Thus:

> f := -8*Pi^2*n*sqrt((1/2)*Pi)*(Int(...));
   evalf(f);

If you already have the Maple worksheet with the output, try Export to..., Maple Text, and enter a file name of your choice.  If you want only output and not input, then in the View menu choose Show/Hide Contents and uncheck Input.

Alternatively, if you want future output sent to a file (and not the worksheet), you can use the writeto command.  Thus:

> writeto("c:/test/myfile.txt");

To return to normal,

> writeto(terminal);

If you use writeto in the Standard interface, make sure Typesetting Level is Maple Standard rather than Extended (use Tools, Options..., Display) ,

AFAIK, the LinearAlgebra package's Eigenvalues and Eigenvectors commands can't do this.  Of course you could take all the eigenvalues and select the ones you want.  Or you could write your own procedure, but it would probably end up being much slower than Eigenvalues.

Because you left out a multiplication sign, your d1 is interpreted as sec(I*(x^2-2)) (which is real-valued). Putting the multiplication sign in, you get an expression that is real only when x is a multiple of Pi.

Now I'm not sure what you mean by "intersection point".  Are you trying to solve d1 = d2?  That won't have any real solutions.  

Another way is to use Sample in the Statistics package.

> with(Statistics):
   U:= RandomVariable(Uniform(0,1));
   S:= sort(Sample(U, 6));
   

In all likelihood there are no closed-form solutions.  Besides series, you could try numerical solutions.  Using Pi, as jakubi remarked:

> S:= dsolve({de, ics}, y(t), numeric);

> plots[odeplot](S,[t,y(t)],t=0..10);

It looks like y(t) ~ 1/(2 t) as t -> infinity.  More precisely:

1/(2*t)+1/(2*Pi*t^3)+3/2/Pi^2/t^5+1/48*(Pi^2+360)/Pi^3/t^7+1/6*(2*Pi^2+315)/Pi^4/t^9+3/1280*(Pi^4+1920*Pi^2+201600)/Pi^5/t^11 + O(t^(-13))

> plot(x^2, x = 0 .. 5, labels = [x, rho*[`kg/m`^3]]);

As for legends in 3D plots: no, that option is not available.  I guess you could simulate a legend to some extent with textplot3d.

It seems the filled option in surfdata is ignored.  You could try something like this.

> with(plots):
   L:= [[1, 1, .69], [1, 2, .48], [2, 1, .37], [2, 2, .44]]; # your data
   m:= nops(L); n:= nops(L[1]);  # number of rows and columns
   toppts:= [seq([seq([i/(m-1),j/(n-1),L[i+1,j+1]],j=0..n-1)],i=0..m-1)]; # points on the top
   basepts:= [seq([seq([i/(m-1),j/(n-1),0],j=0..n-1)],i=0..m-1)];
   perimeter:= [seq([1,j],j=1..n), seq([i,n],i=1+1..m),seq([m,n-j],j=1..n-1),seq([m-i,1],i=1..m-1)];
   sides:= seq(polygonplot3d([toppts[op(perimeter[i])],toppts[op(perimeter[i+1])],
      basepts[op(perimeter[i+1])],basepts[op(perimeter[i])]]),i=1..nops(perimeter)-1):
   base:= polygonplot3d([[0,0,0],[1,0,0],[1,1,0],[0,1,0]]):
   display([surfdata(L), sides, base], axes=frame, labels=[x,y,z]);

There is presumably no closed-form expression for y(r).  If you want a good approximation on a particular interval, you might try something from the numapprox package, perhaps chebyshev or chebpade.   Interpolation is _not_ likely to give good results.  But I don't know why you think you need "an expression" for y(r).

Hmmm, it looks like chebpade has trouble with dsolve(..., numeric) output.  This may warrant further investigation.

> de := diff(y(r), r, r)+2*(diff(y(r),r))/r+9*(16*43)*Pi^2*sqrt(Pi/(2*(1/43)^3))
*polylog(3/2,-exp(1/43*(43-y(r))))/(16*Pi^2*sqrt(43)^3) = 0;

Change of independent variable: s = 1/r.

> sde:= simplify(PDEtools[dchange](r=1/s, de));

sde := s^4*diff(y(s),`$`(s,2))+387/2*2^(1/2)*Pi^(1/2)*polylog(3/2,-exp(1-1/43*y(s))) = 0

As mentioned in another thread, it looks like the solutions are asymptotically logarithmic.  Try a logarithm:

> simplify(eval(lhs(sde),y(s) = c + a*ln(s)));

-s^2*a+387/2*2^(1/2)*Pi^(1/2)*polylog(3/2,-exp(1-c/43)*s^(-1/43*a))

Now polylog(3/2, t) = t + O(t^2), so for this to give 0 to lowest order, we want a = -86.

> series(eval(%,a=-86),s);

(86-387/2*2^(1/2)*Pi^(1/2)*exp(1-1/43*c))*s^2+387/4*Pi^(1/2)*exp(1-1/43*c)^2*s^4+O(s^6)

Choose c to make the s^2 term 0.

> c0:= solve(coeff(%,s,2)=0);

c0 := 43+43/2*ln(81/8*Pi)

Now let y(s) = c0  - 86*ln(s) + u(s), where hopefully u(s) is small when s is small.

> ude:= simplify(eval(sde, y(s) = c0 - 86*ln(s) + u(s)));

ude := 86*s^2+s^4*diff(u(s),`$`(s,2))+387/2*2^(1/2)*Pi^(1/2)*polylog(3/2,-2/9*2^(1/2)/Pi^(1/2)*s^2*exp(-1/43*u(s))) = 0

Linearize this, and solve the linearization:

> lde:= eval(ude,{polylog=unapply(t,(r,t)), exp=(t -> 1+t)});
  dsolve(lde);

u(s) = _C1*s^(1/2)*sin(1/2*7^(1/2)*ln(s))+_C2*s^(1/2)*cos(1/2*7^(1/2)*ln(s))

Presumably the _C1 and _C2 will be arbitrary constants, corresponding to the fact that solutions form a 2-parameter family.  The next terms, which we'll want to make the s^3 terms in the differential equation 0, are, I think, c1*s*cos(7^(1/2)*ln(s))+c2*s*sin(7^(1/2)*ln(s))+c3*s.

> eval(ude,u(s) = _C1*s^(1/2)*sin(1/2*7^(1/2)*ln(s))+_C2*s^(1/2)*cos(1/2*7^(1/2)*ln(s)) 
    + c1*s*cos(7^(1/2)*ln(s)) + c2*s*sin(7^(1/2)*ln(s)) + c3*s);

The usual series command can't handle this, but the MultiSeries version can. 

> MultiSeries:-series(lhs(%),s,3);

(2/43*_C1*sin(1/2*ln(1/s)*7^(1/2))*_C2*cos(1/2*ln(1/s)*7^(1/2))-5*c1*cos(ln(1/s)*7^(1/2))-1/43*_C1^2*sin(1/2*ln(1/s)*7^(1/2))^2-1/43*_C2^2*cos(1/2*ln(1/s)*7^(1/2))^2+5*c2*sin(ln(1/s)*7^(1/2))+7^(1/2)*c2*cos(ln(1/s)*7^(1/2))+2*c3+7^(1/2)*c1*sin(ln(1/s)*7^(1/2)))*s^3+O(s^(7/2)*(((-1/86*exp(-I*ln(1/s)*7^(1/2))-1/86*exp(ln(1/s)*7^(1/2)*I))*c1+(1/86*I*exp(-I*ln(1/s)*7^(1/2))-1/86*I*exp(ln(1/s)*7^(1/2)*I))*c2-1/43*c3)*((1/86*I*exp(-1/2*I*ln(1/s)*7^(1/2))-1/86*I*exp(1/2*I*ln(1/s)*7^(1/2)))*_C1+(-1/86*exp(-1/2*I*ln(1/s)*7^(1/2))-1/86*exp(1/2*I*ln(1/s)*7^(1/2)))*_C2)+1/6*((1/86*I*exp(-1/2*I*ln(1/s)*7^(1/2))-1/86*I*exp(1/2*I*ln(1/s)*7^(1/2)))*_C1+(-1/86*exp(-1/2*I*ln(1/s)*7^(1/2))-1/86*exp(1/2*I*ln(1/s)*7^(1/2)))*_C2)^3))

> s3term:= simplify(combine(eval(%,O=0))) assuming s > 0;

s3term := 1/86*s^3*(-2*_C1*_C2*sin(7^(1/2)*ln(s))-430*c1*cos(7^(1/2)*ln(s))-_C1^2+_C1^2*cos(7^(1/2)*ln(s))-_C2^2*cos(7^(1/2)*ln(s))-_C2^2-430*c2*sin(7^(1/2)*ln(s))+86*7^(1/2)*c2*cos(7^(1/2)*ln(s))+172*c3-86*7^(1/2)*c1*sin(7^(1/2)*ln(s)))

> solve(identity(s3term,s),{c1,c2,c3});

{c1 = 1/19264*7^(1/2)*(-5*7^(1/2)*_C2^2-14*_C1*_C2+5*7^(1/2)*_C1^2), c2 = -1/19264*(-7*_C2^2+10*_C1*_C2*7^(1/2)+7*_C1^2)*7^(1/2), c3 = 1/172*_C2^2+1/172*_C1^2}

To sum up, an asymptotic form for the solution as r -> infinity is

with c0, c1, c2 and c3 as above.  The approximate values of _C1 and _C2 may be obtained by comparing this function and its derivative to the values for the numerical solution (obtained in another thread) at some convenient large value of r.

You might try this:

> seqindets:= proc(X,T)
   if type(X,T) then X
   elif type(X,atomic) then NULL
   else seq(seqindets(x,T),x=X)
   fi
 end proc;

> [seqindets]( [[`@`(1, {.4}), `@`(2, {.4})], [[], 1], [[], 2], 3],numeric);

[1, .4, 2, .4, 1, 2, 3]

You may get better results by solving a larger system of differential equations.  Thus to your system involving c(t), add the differential equation

diff(U(t),t) = exp(-rho*t)*c(t)^(1-theta)/(1-theta)

and boundary condition U(0) = 0.  Your Uzero will then be U(400).

The "9 months after a holiday" idea is cute, but not really supported by statistics I think.  At <http://www.cdc.gov/nchs/data/statab/natfinal2003.annvol1_16.pdf> you can find some actual statistics: the number of live births in the US for each day of 2003.  The most obvious pattern is that there are relatively few births on weekends and holidays, presumably because doctors would not schedule an elective caesarian or induction for those days.  There are also seasonal effects: The least popular birthday is December 25 (6628 births), the most popular is December 30 (14438).