Your first error is that the assignment operator in Maple is ":=", not "=".

You seem to want a system of two differential equations in z_1 and z_2, but you're only giving dsolve one equation at a time.

dsolve can solve complex systems numerically, but you must tell it about the derivative of z_i, not the derivative of conjugate(z_i).  So you should try something like

> dsolve({diff(z_1(t),t) - conjugate(summ_1) = 0, diff(z_2(t),t) - conjugate(summ_2) = 0,
      z_1(0) = I, z_2(0) = 2*I}, numeric);

For example, you might try:

for n to 20 do 'A'^n = A^n end do;

Assuming your "natural numbers" don't include 0, the generating function for the number of solutions of 2x+3y+5z+8t=n is

> convert(s^2/(1-s^2)*s^3/(1-s^3)*s^5/(1-s^5)*s^8/(1-s^8),FPS,s);

So you want to evaluate

at k = 10^2012.  Of course the result is rather big.  The roots of Z^4 + Z^3 + Z^2 + Z + 1 being 5th roots of 1, the sum over those roots is periodic with period 5; that value for k = 0 mod 5 is -1/5.  Similarly, the sum over the roots of Z^2 + 1 when k = 0 mod 4 is -1/16, the sum over the roots of Z^4 + 1 when k = 0 mod 8 is -1/4, and the sum over the roots of Z^2 + Z + 1 when k = 1 mod 3 is -1/9.  Evaluating the rest at k = 10^2012 is easy.  The final result can be seen in the following worksheet.  It is an integer of 6033 digits starting with 6944 and ending with 1110.

drops.mw

For a dense Matrix containing floats, if hardware precision is ok, it uses NAG routies F08KEF and F08MEF.  F08KEF reduces the Matrix to bidiagonal form using orthogonal matrices.  If only the singular values are required, F08MEF computes them using the differential qd algorithm.  See
www.nag.com/numeric/Fl/manual/pdf/F08/f08kef.pdf and
www.nag.com/numeric/Fl/manual/pdf/F08/f08mef.pdf.

Instead of Mean(exp(a*X)), which may involve a difficult sum or integral, it's probably better to use MGF(X,a), especially in cases where the distribution of X is "well-known".  For example, on my computer:

> with(Statistics):
assume(a>0, a 1, t < 1); X:= RandomVariable(GammaDistribution(a,b));
time(Mean(exp(X*t)));

  2.324

> restart; with(Statistics):
assume(a>0, a 1, t < 1); X:= RandomVariable(GammaDistribution(a,b));
time(MGF(X,t));

0.046

Basically, if you're interested in speed and memory usage, you want to avoid as much as possible doing difficult symbolic computations on-the-fly.  In this case, I expect that MGF is just doing a look-up and finding a known formula.

Your ODE's are almost OK, but your initial conditions are not.  You have

I don't know how y(1) gets into the first equation: I suspect it should be y1(t).

Initial conditions should be at t=0, not at arbitrary t.  If y1(0), y2(0) and y3(0) are given, that would determine the second derivatives at t=0.  The initial conditions should instead specify the first derivatives at t=0, thus D(y1)(0) = ..., D(y2)(0) = ..., D(y3)(0) = ....

Try this:

> RNN:= proc(M::{Array,Matrix},k::posint, t::type :=numeric)
  local U,d,i,dk;
  d:= [rtable_dims(M)];
  U:= remove(r -> hastype(M[op(subsop(k=r,d))],Non(t),'exclude_container'),[seq(i,i=d[k])]);
  M[op(subsop(k=U,d))]
end proc;
     

This will remove the part of M where there is a non-numeric entry sharing the same k'th index. Thus for a Matrix or 2-dimensional Array, RNN(M, 1) returns the rows that are all numeric, RNN(M, 2) returns the columns that are all numeric.  

An optional third argument replaces numeric by any other type.  Thus RNN(M, 1, positive) returns rows where all entries are positive.

> with(VectorCalculus):
   SetCoordinates(cartesian[x,y,z]);
   V:= VectorField(<u(x,y,z),v(x,y,z),w(x,y,z)>);
   (V . Nabla)~(V);

This should work in Maple 14 or 15.  In previous versions, before the elementwise operator ~ was introduced, you could use

> map(V . Nabla, V);

While Markiyan is right that simplify(..., symbolic) works, that shouldn't be necessary here: this simplification should be valid for all complex numbers x (except, of course, the square roots of -3/2).  If it were not always valid, appropriate assumptions should work; but in this case I can't get anything of the form
simplify(..) assuming ... to work.  Note that if the sqrt(4*x^2+6) was replaced by sqrt(2)*sqrt(2*x^2+3), the expression would automatically be simplified to

and simplify would do the rest.

I'll submit an SCR.

In Cartesian coordinates (for an incompressible Newtonian fluid):

> with(VectorCalculus):
   SetCoordinates(cartesian[x,y,z]):
   V:= VectorField(< seq(v[i](x,y,z,t),i=1..3)>);
   NavierStokes:= rho*(diff(V,t) + (V . Nabla)~(V)) =        
        - Gradient(p(x,y,z,t)) + mu*Laplacian~(V) + VectorField(< seq(f[i](x,y,z,t),i=1..3)>);

Similarly in conical coordinates, but you'd use

   SetCoordinates(conical[u,v,w]);

and replace x,y,z by u,v,w.

For example, if your data file "c:/text.txt" looks like this:

1      3.6
2.3   3.9        
4.5   4.3
6.7   4.5
8.9   4.6

(where each line consists of an x value and the corresponding y value, separated by spaces)
and you want to fit this to

you could do it this way.

> data:= readdata("c:/text.txt",2);
   Statistics:-Fit(exp(a+b*x+c*x^2),data,x);

You can do a 3D plot of your solution directly with e.g.

> pds:-plot3d(t=0..1);

As for plotting an Array, you could do that with

> plots[surfdata](Array(A, datatype=float[8]), axes=box);

The datatype=float[8] option could have been included in the original definition of A, so it would be

> A := Array([seq([seq(U((1/20)*j, (1/20)*k), j = 0 .. 20)], k = 0 .. 20)],datatype=float[8]);
  plots[surfdata](A, axes=box);

A just contains the values U(j/20, k/20) for j from 0 to 20 and k from 0 to 20, so A[1,1] = U(0,0).  If you wanted U(-10, 0), you could have taken j from -20 to 20.

For example:

> eqs:= [ 3*diff(x(t),t$2) + 4*diff(y(t),t) + 2*y(t) + 2= 0,
    5*diff(y(t),t$2)+6*diff(x(t),t) + 7*x(t) + 3*sin(t)= 0];
   vars:= [x(t),y(t)];
   dvars:= [ op(diff~(vars,t$2)), op(diff~(vars,t)),op(vars)];
   M,d:= LinearAlgebra[GenerateMatrix](eqs,dvars);

> A,B,C:= M(.., 1..2), M(..,3..4), M(..,5..6);

The current version is Maple 15.  See http://maplesoft.com/products/Maple/index.aspx

You might look at densityplot in the plots package.  Thus:

> plots[densityplot](f,-4..4, -4..4, grid=[50,50],style=patchnogrid);