Boundary conditions for what?  What are you going to do with them?

If you're dealing with a finite number of functions, say F[m,n] for m from 0 to M and n from 0 to N, you could try something like

seq(seq(F[m,n](1) = `if``(m=0 and n=0, 1,  0), m=0..M), n=0..N)

(and similarly for F[m,n](-1)...)

> solve(2*x-a/x+a-x+3-x=1-(x^2)+3, x);

For any continuously differentiable function f, the differential equation

dy/dx = 1/f'(y)

has solutions f(y) = x + C, i.e. y = f^(-1)(x + C) for arbitrary constant C.

One important difference between a PDE and an ODE is that while the solutions of an ODE will involve arbitrary constants, the solutions of a PDE will typically involve arbitrary functions.  In this case Maple is telling you that there are solutions that are arbitrary  functions of x^2 + y^2.  You could replace _F1 with any (sufficiently differentiable) function of one variable.

Use the RowOperation command in the LinearAlgebra package.

I hope I understand what you're asking for.

> for i from 1 to 2 do
       e[i] := f[i](x[1],x[2],y[1],y[2]);
       X[i]:= x[i](y[1],y[2]);
   end do;
   S:= solve({seq(seq(diff(subs(x[1]=X[1],x[2]=X[2], e[i]), y[j]), i=1..2), j=1..2)},
       {seq(seq(diff(X[i],y[j]), i=1..2), j=1..2)});
   for i from 3 to 4 do
     e[i]:= subs(S, X[1]=x[1], X[2]=x[2], diff(g[i-2](X[1],X[2],y[1],y[2]), y[i-2]))
   end do;

And then solve the system of four equations e[1]=0, e[2] = 0, e[3]=0, e[4] = 0.

You have the first part of the trajectory, but it leaves the box you specified, with y(t) = 20 at approximately t = 18.55.

It looks to me like r(t) -> infinity, x(t) -> -infinity, y(t) -> infinity and z(t) -> 0 as t -> infinity, so you're not heading to any attractor.  Note that with z(t) = 0 the equations for x,y,r have a solution r(t) = exp(t/20), x(t) = -100/79 * exp(t/20), y(t) = 105/79*exp(t/20),
and when x(t) -> -infinity the equation for z(t) will drive z(t) to 0.

.

On the other hand, it's not 0 in general.

> simplify(collect(eq,s,factor)) assuming s>1;

       0

> plot3d(abs(eq),s=-2..2,t=-2..2,axes=box,shading=zhue);

The way to deal with a plot (parametric or not) that "blows up" is to use the view option.  For example:

> plot([t + 1/t, t - 1/t, t = -5 .. 5], view=[-5..5, -5 .. 5], discont=true);

I'll look at the planar case.

The three parameters only occur in two combinations Dd/alphaa and Dcj/alphaa, so I'll write Dd/alphaa = b and Dcj/alphaa = c.  Moreover, the dependent variable occurs only as derivatives, so write diff(zs(r),r) = v(r).  The differential equation becomes

> diff(v(r),r) = (1+v(r)^2)*b-(1+v(r)^2)^(3/2)*c;

This separable differential equation has the implicit solution

> r = int(1/(1+v^2)/(b-c*sqrt(1+v^2)), v) + C;

The integral can be done in closed form, but it's complicated.  Assuming b and c are positive:

> simplify(%) assuming b>0,c>0;

r = 1/2*1/b/(b^2-c^2)^(1/2)*(c*ln(-(v*(b^2-c^2)^(1/2)+c+b*(1+v^2)^(1/2))/(-c*v+(b^2-c^2)^(1/2)))-c*ln((-v*(b^2-c^2)^(1/2)+c+b*(1+v^2)^(1/2))/(c*v+(b^2-c^2)^(1/2)))+2*arctanh(c*v/(b^2-c^2)^(1/2))*c+2*(b^2-c^2)^(1/2)*arctan(v)+2*(b^2-c^2)^(1/2)*C*b)

It seems highly unlikely that this can be solved explicitly for v in terms of r.  For example, with b = c = 1:

> r = int(1/(1+v^2)/(1-sqrt(1+v^2)), v) + C;

r = 1/v+arctan(v)+1/v*(1+v^2)^(1/2)+C

I think the sum does not have a closed-form expression. 

If you write exp(-c/x) = t, and assume c/x > 0 (so 0 < t < 1) and |z| < 1, your sum becomes

Sum(z^j/(1-t^j), j=1..infinity)

which can be written as a double sum

Sum(Sum((z*t^k)^j, k=0..infinity), j=1..infinity)

Interchanging the order of summation, this is

Sum(z*t^k/(1-z*t^k), k=0..infinity)

but that also offers little help.

You could get numerical solutions for given values of t and N, or series solutions in powers of N or in powers of t.

Is this the sort of thing you're interested in?

> p1:= 1+Sum(c[i]*e^i, i=1..infinity);
  p2:= 1 + p1;
  map(normal,series(p1/p2, e, 7));
1/2+1/4*c[1]*e+(1/4*c[2]-1/8*c[1]^2)*e^2+(1/4*c[3]-1/4*c[1]*c[2]+1/16*c[1]^3)*e^3+(-1/8*c[2]^2+3/16*c[2]*c[1]^2-1/4*c[1]*c[3]+1/4*c[4]-1/32*c[1]^4)*e^4+(-1/4*c[3]*c[2]+3/16*c[3]*c[1]^2-1/4*c[1]*c[4]+1/4*c[5]+3/16*c[1]*c[2]^2-1/8*c[2]*c[1]^3+1/64*c[1]^5)*e^5+(-1/4*c[4]*c[2]+3/16*c[4]*c[1]^2+1/4*c[6]-1/4*c[1]*c[5]-1/8*c[3]^2+3/8*c[3]*c[1]*c[2]-1/8*c[3]*c[1]^3+1/16*c[2]^3-3/16*c[2]^2*c[1]^2+5/64*c[2]*c[1]^4-1/128*c[1]^6)*e^6+O(e^7)

You can create a .wav file using the AudioTools package.  Let's say you have a list L of notes to play, and each is to have a duration of 0.5 seconds.  For example,

> L := [3,2,1,2,3,3,3,2,2,2,3,5,5]; # Mary had a little lamb

> with(AudioTools):
   n:= nops(L);
   A:= Create(duration=0.5*n); 
   jmax := op(2, ArrayDims(A));

Now with middle C corresponding to a frequency of  about 261.6 Hertz, and each semitone above that multiplies the frequency by 2^(1/12).  So here are the frequencies for notes 0 to 7 (note that E to F is a semitone, while the other intervals are two semitones:

> freq[0]:= 261.6:
  for i from 1 to 7 do
      if i = 4 then freq[i] := freq[i-1]*2^(1.0/12)
      else freq[i]:= freq[i-1]*2^(2.0/12)
      fi
   od:

I'll use a pure sine-wave.  The signal corresponding to note k is thus f(t) = sin(2*Pi*freq[k]*t).  Each time step corresponds to a time of 0.5*n/jmax seconds.  Time step #j is note L[i] where i = ceil(j*n/jmax).  To avoid notes running into each other, I'll leave a  small gap between them

> for j from 1 to jmax do
      if frac(j/jmax*n) < 1/10 then A[j] = 0 else
      A[j]:= evalf(sin(2*Pi*freq[L[ceil(j*n/jmax)]]*j*0.5*n/jmax))
   fi od;

Now to write the file:

> Write("c:/temp/mywav.wav", A);

Not very musical, but it works.

Let B(x) = B0+B1*x+B2*x^2+B3*x^3+B4*x^4+B5*x^5+B6*x^6+B7*x^7,
A(x) = a0+a1*x+a2*x^2+a3*x^3+a4*x^4 and C(x) = c0+c1*x+c2*x^2+c3*x^3.

So you want B(x) * C(x) to be a power of A(x).  Assuming B7 and c3 are nonzero,  B(x)*C(x) would be a polynomial of degree 10.  On the other hand, if a4 is nonzero A(x) has degree 4,
so powerN must be 10/4 = 5/2.  Then A(x) must be the square of a quadratic polynomial.  Now this quadratic has one root of multiplicity 2 or two roots of multiplicity 1, making B(x)*C(x)  have either one root of multiplicity 10 or two roots of multiplicity 5 each: 7 of these 10 belong to B(x) and the other 3 to C(x).  Factoring B(x) should tell what the situation is.

For example:

> B := x^7-9*x^6+9*x^5+135*x^4-405*x^3-243*x^2+2187*x-2187;
   factor(B);

(x-3)^5*(x+3)^2

Thus we have the two distinct roots -3 and +3.  We must have (up to a constant)
C(x) = (x+3)^3 and A(x) = (x-3)^2 * (x+3)^2.

I'll leave the other cases up to you.

Note that your system of DE's is second order in y and z, but only first order in x.  So you may have too many boundary conditions.  It's also not at all clear that you can put on those conditions at infinity, even if they would make sense at a finite t.