I don't understand why pdsolve gives this "incorrect number of boundary conditions" error with the new equation.  The same derivatives occur in both the old Aeq and the new Aeq.  Maybe someone more expert in PDE's can enlighten us, or maybe it's a bug. For reference, the old Aeq was

diff(diff(diff(a(t,y),t),y),y) = (-2*a(t,y)*diff(a(t,y),t)*diff(diff(a(t,y),y),y)^2+2*diff(a(t,y),t)*diff(diff(a(t,y),y),y)
*diff(a(t,y),y)^2-2*diff(a(t,y),t)*diff(a(t,y),y)^4*a(t,y)+2*diff(a(t,y),t)*diff(a(t,y),y)^4-6*diff(a(t,y),t)
*diff(diff(a(t,y),y),y)*diff(a(t,y),y)^2*a(t,y)-a(t,y)*diff(a(t,y),t)*diff(a(t,y),y)*diff(diff(diff(a(t,y),y),y),y))
/(diff(a(t,y),y)^2*a(t,y)+a(t,y)*diff(diff(a(t,y),y),y))

and the new one is

diff(diff(diff(a(t,y),t),y),y) = (-3*a(t,y)*diff(a(t,y),t)*diff(diff(a(t,y),y),y)^2-a(t,y)*diff(a(t,y),t)*diff(a(t,y),y)
*diff(diff(diff(a(t,y),y),y),y)-6*diff(a(t,y),t)*diff(diff(a(t,y),y),y)*diff(a(t,y),y)^2)/(a(t,y)^2*diff(diff(a(t,y),y),y)
+diff(a(t,y),y)^2*a(t,y))

and in both cases the boundary and initial conditions are

{a(0,y) = y, a(t,0) = 0, a(t,1) = 1+2*t, D[2](a)(t,0) = 1+t}

I don't understand why pdsolve gives this "incorrect number of boundary conditions" error with the new equation.  The same derivatives occur in both the old Aeq and the new Aeq.  Maybe someone more expert in PDE's can enlighten us, or maybe it's a bug. For reference, the old Aeq was

diff(diff(diff(a(t,y),t),y),y) = (-2*a(t,y)*diff(a(t,y),t)*diff(diff(a(t,y),y),y)^2+2*diff(a(t,y),t)*diff(diff(a(t,y),y),y)
*diff(a(t,y),y)^2-2*diff(a(t,y),t)*diff(a(t,y),y)^4*a(t,y)+2*diff(a(t,y),t)*diff(a(t,y),y)^4-6*diff(a(t,y),t)
*diff(diff(a(t,y),y),y)*diff(a(t,y),y)^2*a(t,y)-a(t,y)*diff(a(t,y),t)*diff(a(t,y),y)*diff(diff(diff(a(t,y),y),y),y))
/(diff(a(t,y),y)^2*a(t,y)+a(t,y)*diff(diff(a(t,y),y),y))

and the new one is

diff(diff(diff(a(t,y),t),y),y) = (-3*a(t,y)*diff(a(t,y),t)*diff(diff(a(t,y),y),y)^2-a(t,y)*diff(a(t,y),t)*diff(a(t,y),y)
*diff(diff(diff(a(t,y),y),y),y)-6*diff(a(t,y),t)*diff(diff(a(t,y),y),y)*diff(a(t,y),y)^2)/(a(t,y)^2*diff(diff(a(t,y),y),y)
+diff(a(t,y),y)^2*a(t,y))

and in both cases the boundary and initial conditions are

{a(0,y) = y, a(t,0) = 0, a(t,1) = 1+2*t, D[2](a)(t,0) = 1+t}

I don't understand what you want to do with max.  Do you mean you want the member(s) of the list a where the second entry is maximized?

> a:=[ [ 6,3],[5,6],[6,6],[5,8],[4,7],[5.5,10],[6,12],[5,11],[5,13]]:
  m:= max(map2(op,2,a)):  
  select(t -> (t[2] = m), a);

I don't understand what you want to do with max.  Do you mean you want the member(s) of the list a where the second entry is maximized?

> a:=[ [ 6,3],[5,6],[6,6],[5,8],[4,7],[5.5,10],[6,12],[5,11],[5,13]]:
  m:= max(map2(op,2,a)):  
  select(t -> (t[2] = m), a);

No, it doesn't: a, b and c are assigned the same Array, so that a change in one affects the others.  The same happens using map:

> L:= [a,b,c]:  
  assign(map(`=`,L,Array(1..16)));
  a[1]:= foo:
  b[1];

foo

The solution with ~ seems OK though, because you're making three distinct Arrays.

No, it doesn't: a, b and c are assigned the same Array, so that a change in one affects the others.  The same happens using map:

> L:= [a,b,c]:  
  assign(map(`=`,L,Array(1..16)));
  a[1]:= foo:
  b[1];

foo

The solution with ~ seems OK though, because you're making three distinct Arrays.

> A:= a(t,y): N:= n(t,y):
  sys:= {diff(A,t)^2 - N^2*diff(A,y$2) - N^2*diff(A,y)^2 = 0,
   N^3*diff(A,y)^2 + N^2*A*diff(A,y)*diff(N,y) -N*A*diff(A,t$2)- N*A*diff(A,t)^2+A*diff(A,t)*diff(N,t)=0};
  C:= PDEtools[casesplit](sys);
  Aeq:= op([1,2],C[1]);
  ibc2:= {a(0,y)=y,a(t,0)=0,D[2](a)(t,0)=1+t,a(t,1)=1+2*t};
  Sol2:= pdsolve(Aeq,ibc2,numeric);
  af:= subs(Sol2:-value(output=listprocedure),a(t,y));
  nf:= D[1](af)/sqrt(D[2,2](af)+D[2](af)^2);
  plot3d(nf, 0..1, 0..1, axes=box);

I think this should work (but it looks like it will take a lot of time).

> A:= a(t,y): N:= n(t,y):
  sys:= {diff(A,t)^2 - N^2*diff(A,y$2) - N^2*diff(A,y)^2 = 0,
   N^3*diff(A,y)^2 + N^2*A*diff(A,y)*diff(N,y) -N*A*diff(A,t$2)- N*A*diff(A,t)^2+A*diff(A,t)*diff(N,t)=0};
  C:= PDEtools[casesplit](sys);
  Aeq:= op([1,2],C[1]);
  ibc2:= {a(0,y)=y,a(t,0)=0,D[2](a)(t,0)=1+t,a(t,1)=1+2*t};
  Sol2:= pdsolve(Aeq,ibc2,numeric);
  af:= subs(Sol2:-value(output=listprocedure),a(t,y));
  nf:= D[1](af)/sqrt(D[2,2](af)+D[2](af)^2);
  plot3d(nf, 0..1, 0..1, axes=box);

I think this should work (but it looks like it will take a lot of time).

Here is an animation.  The red, blue and green curves intersect at solutions.  It looks like for J > 4 there is one solution inside the box -1 < x,y,z < 1.

> with(plots): 
  Frame:= proc(j) 
    local R; 
    R:= eval([R1,R2,R3],J=j);
    display([intersectplot(R[1],R[2],x=-1..1,y=-1..1,z=-1..1,colour=red),
             intersectplot(R[1],R[3],x=-1..1,y=-1..1,z=-1..1,colour=green),
             intersectplot(R[2],R[3],x=-1..1,y=-1..1,z=-1..1,colour=blue)])
  end proc; 
  animate(Frame, [j], j=1..10, axes=box);

 (This will take a while, and it's late - I'll post a link to an animated .gif tomorrow)

The following three polynomials must be 0 for a solution:

[256-192*y^2*J^2*z^2*x^2-1024*y-16*y^4*J^3*z^3+96*y^4*J^2*z^2-256*y^4*J*z+16*y^3*J^3*z^3-192*y^3*J^2*z^2+768*y^3*J*z+320*y^3*x^2*J^2*z^2-1280*y^3*x^2*J*z-16*y^3*J^3*z^3*x^2+16*y^4*J^3*z^3*x^2-128*y^4*J^2*z^2*x^2+128*y^2*J^2*z^2*x^4+256*y^4*J*z*x^2+64*y^3*J^2*x^2+512*y*J*z*x^4-1024*y^2*J*z*x^4+512*y^3*x^4*J*z-128*y^3*x^4*J^2*z^2-2*y^4*J^4*x^2*z^2+16*y^4*J^3*x^2*z-16*y^3*J^3*x^2*z+y^4*J^4*z^4+16*y^3*J^3*x^4*z-16*y^4*J^3*x^4*z-3072*x^2*y^2+3072*y*x^2+96*y^2*J^2*z^2-768*y^2*J*z+256*y*J*z-1024*y^3-1024*x^2+1536*y^2+256*y^4-32*y^2*J^2*x^2+1792*y^2*J*z*x^2-768*y*J*z*x^2+1024*y^2*x^4+1024*x^4+1024*y^3*x^2-2048*y*x^4+y^4*x^4*J^4-32*y^4*J^2*x^2+64*y^4*x^4*J^2-64*y^3*x^4*J^2, 64*J^2*x^2*z^2-512*y^2*J*x-320*y^2*J^2*z^2*x^2-128*y^4*J^2*z^2+512*y^4*J*z+192*x^3*z*y^4*J^2-16*z^2*J^3*x^3*y^4-192*J^2*y^4*x*z+16*J^3*z^2*x*y^4-16*J^3*z^4*x*y^4+192*J^2*y^4*x*z^3-4*J^4*y^2*x^2*z^2+32*x^2*J^3*z^3*y^2-256*J^2*y^2*x*z^3+32*J^3*z^4*x*y^2+J^4*z^4+32*J^2*z^4+16*x^2*J^3*z+16*J^3*z^2*x-64*J^2*x*z-2*J^4*y^4*z^2+64*y^4*J^2*z^4-512*y^4*J*z^3+64*x^3*z*J^2+512*y^2*J*z^3-96*y^2*J^2*z^4-16*z^2*J^3*x^3+2*J^4*x^2*z^2-2*J^4*y^2*z^4+64*J^2*x*z^3-16*J^3*z^4*x-16*x^2*J^3*z^3-16*x^4*J^3*z-32*J^2*z^2-2*J^4*z^2+4*J^4*y^2*z^2-16*y^4*J^3*z^3*x^2+256*y^4*J^2*z^2*x^2-64*y^2*J^2*z^2*x^4-1024*y^4*J*z*x^2+256*y^2*J*z*x^4-2*y^4*J^4*x^2+1024*z^2*J*x*y^2+2*y^4*J^4*x^2*z^2+16*y^4*J^3*x^2*z-128*x^3*J^2*z^3*y^2+768*x^3*J*z^2*y^2+768*x^2*J*z^3*y^2-512*y^4*J*x^3+512*x^3*J*y^2-64*J^2*x^2*z^4*y^2-1024*y^4*z^2*J*x+256*z^4*J*x*y^2-32*J^3*z^2*x*y^2+256*J^2*y^2*x*z+32*z^2*J^3*x^3*y^2-256*x^3*z*y^2*J^2+32*x^4*J^3*z*y^2-32*x^2*J^3*z*y^2+y^4*J^4*z^4+64*J^2*x^2*z^4-256*x^4*J*z+128*x^3*J^2*z^3-768*x^3*J*z^2-768*x^2*J*z^3+64*x^4*J^2*z^2-2*J^4*x^2-256*z^4*J*x-16*y^4*J^3*x^4*z+4*J^4*y^2*x^2-2*x^4*J^4*y^2+512*y^4*J*x-64*y^2*J^2-1024*x^2*y^2+160*y^2*J^2*z^2-512*y^2*J*z+1024*x^3*z+1024*x*z^3-1024*y^2*z^2-1024*y^2*x^3*z+256*z^4+160*y^2*J^2*x^2+1024*y^2*J*z*x^2+1024*y^4*x^2+2048*y^4*x*z+1024*y^4*z^2-1024*y^2*x*z^3-32*x^2*J^2+32*x^4*J^2+x^4*J^4-2048*y^2*x^2*z^2+1536*x^2*z^2+256*x^4+J^4-2*J^4*y^2+y^4*J^4-2048*y^2*x*z+64*y^4*J^2+y^4*x^4*J^4-128*y^4*J^2*x^2+64*y^4*x^4*J^2-96*x^4*y^2*J^2, 256+512*y^2*J*z*x-256*y*J*z*x+64*y^6*J^2*x^2+128*y^2*J^2*z^2*x^2-1024*y+256*y^5*z*J*x-64*y^5*J^2*z^2-512*y^4*J*z*x-64*y^3*x^2*J^2*z^2-192*y^4*J^2*z^2*x^2+64*y^3*J^2*x^2-16*y^5*J^3*z*x+16*y^7*J^3*z*x+16*y^3*J^3*z*x^3+512*z^3*y*J*x-16*y^5*J^3*z*x^3-512*y^3*J*z^3*x+64*y^6*J^2*z^2-128*y^5*J^2*x^2+128*y^5*J^2*x^2*z^2+512*y^4*J*z^3*x-512*y^2*J*z^3*x+y^8*J^4-1024*y^3-3072*y^2*z^2+64*y^5*J^2-32*y^6*J^2+1536*y^2-1024*z^2+256*y^4+1024*z^4-2*x^2*y^6*J^4-32*y^2*J^2*x^2+3072*y*z^2+1024*y^2*z^4+1024*y^3*z^2-2048*y*z^4-32*y^4*J^2+y^4*x^4*J^4+32*y^4*J^2*x^2]

In principle it is possible to eliminate two of the variables, say y and z, leaving one big polynomial in x and J.  But such a polynomial is likely to be enormous.  So "closed-form" solutions are not a realistic hope.  

Of course, for particular values of J you can find solutions numerically.

Another solution is x = -1, y = z = 4/J. 

However (at least for some values of J) there are other solutions, which seem more difficult to find.  For example, at J = 4 there is a solution approximately x = -.9125038849, y = .7622519735, z = .9603396146.

You're missing a multiplication sign after the first y in R2, but even after correcting that your system may take a long time to solve.  I don't know if there's a way to fix this: nonlinear algebraic systems can be hard to solve symbolically, and the solutions can be extremely complicated.

If a, m or n is a float or a fraction, || will not work.  However, you can use sprintf instead.  Thus:

plotsetup(ps, plotoutput=sprintf("c:/tmp/A[%a]M[%a]N[%a].eps", a, m, n), ...

If a, m or n is a float or a fraction, || will not work.  However, you can use sprintf instead.  Thus:

plotsetup(ps, plotoutput=sprintf("c:/tmp/A[%a]M[%a]N[%a].eps", a, m, n), ...

There are links, but they don't seem to work:  I just get "Page Not Found".