Another method for this PDE system would be power series.  Assuming analyticity in a neighbourhood of the origin, we can write

u := sum(sum(a[i,j]*x^i*t^j,j = 0 .. infinity),i = 0 .. infinity)

v := sum(sum(j*a[i,j]*x^i*t^(j-1),j = 0 .. infinity),i = 0 .. infinity)

and then the second PDE says

 a[i,j] = (i+1)*(j+1)*a[i+1,j+1]

so that if m = min(i,j),

a[i,j]=a[i-m,j-m]*(i+j-2*m)!/i!/j!

(that should say a[i,j] = a[i-m,j-m]*(i+j-2*m)!/i!/j! )

Summing over the diagonals

u = a[0,0]*BesselI(0,2*x^(1/2)*t^(1/2))+sum(BesselI(j,2*x^(1/2)*t^(1/2))*j!*(a[0,j]*t^(j/2)/x^(j/2)+a[j,0]*x^(j/2)/t^(j/2)),j = 1 .. infinity)

where

u(x,0) = sum(a[i,0]*x^i,i = 0 .. infinity)

u(0,t) = sum(a[0,j]*t^j,j = 0 .. infinity)

and thus

v(0,t) = sum(a[0,j]*j*t^(j-1),j = 1 .. infinity)

1) Square brackets [] are never used for grouping: they are for lists and indices.  Use parentheses () instead.

2) You have initial conditions but no boundary conditions.  You need two boundary conditions.  Try e.g.:

> IBC = {u(x, 0) = tanh(x), v(x, 0) = 0, u(0,t) = 0, u(1,t) = tanh(1-t)};

I don't think the geom3d package is useful here.
The simplest way to define a 3d arc is using a parametric curve, giving the coordinates as a function of a parameter t in some interval.  For example:

A := <cos(t), sin(t), 0>;    

Rotations can then be done using, e.g.,  RotationMatrix in the Student[LinearAlgebra] package.

Since there's no derivative wrt x, this PDE system is really an ODE system in t with a parameter x. 

diff(u(t),t)=diff(v(t),t,t) + u(t) - 1/3*u(t)^3 - v(t)

diff(v(t),t)=u(t)-v(t)

This turns out to be only first order: substitute u(t) = diff(v(t),t)+v(t) into the first equation and simplify, and you get 1/3*(diff(v(t),t)+v(t))^3 = 0

So you should have only one initial condition for each x, and you have two.

For example:

> with(Maplets[Elements]):
  A:= x^3; B:= x^3;
myProc := proc()
   global A,B;
   A:= int(A,x); B:= diff(B,x);
   Maplets:-Tools:-Set('MV1' = A );
   Maplets:-Tools:-Set('MV2' = B );
end proc:

maplet1 := Maplet( 
[[MathMLViewer[MV1]('value'=A), 
  MathMLViewer[MV2]('value'=B)],
 [Button("Go", Evaluate('function' = "myProc")),
  Button("Done", Shutdown())
 ]]):
Maplets[Display](maplet1);

I don't think your system is the type that the numerical PDE solver can handle.  One independent variable must be a "time" variable, and it must be possible to interpret the system as describing the evolution of the dependent variables in time.  I'm not sure what the precise technical requirements are, but e.g. you have one equation that has no derivative wrt t and another that has no derivative wrt x.

Perhaps if you could post the actual equations we might have some more concrete suggestions on whether other solutions are likely, or other ways to approach the problem. 

You could use the piecewise function, like this:

> R := piecewise(a=a1, b1, a=a2, b2, a=a3, b3);

Your procedure RK2 has x and t as local variables.  Now one line of the procedure is

F1 := unapply(fct, [x, t]);

The x and t in this command will be the local x and t, not the global ones which your fct contains.  You might use :-x and :-t to refer to the global x and t, but I think it would be simpler to  just use different names  for the local variables (maybe X and T).

You don't really have "sets" in the Maple sense, since the orders seem to be significant: in Maple terminology  Y and B could be lists or tables or Arrays.  Rather than embedding the index in the name (e.g. Y1), we indicate the index with square brackets (e.g. Y[1]).  On the other hand, square brackets are _not_ used for grouping: that's the job of parentheses.  

You might also note that (with Q = ln((Y[j]+1)/Y[j]) > 0)

Q*(x-X[k]) >= 1 if and only if x >= X[k] + 1/Q

The basic idea of a numerical solution is that you go from the initial value (in this case at x=0) to whatever value of the independent variable x you're interested in, following what happens to the solution at each stage.  If the solution goes off to infinity as x approaches a certain value, then the numerical solution process can't go past that.  One thing you could try, however, is a transformation.  If z(x) = 1/y(x), then a pole of y(x) corresponds to a zero of z(x).  The numerical solution of the differential equation for z(x) (after cancelling out a factor) won't mind at all having z pass through 0.

The differential equation for z(x):

> zde:= numer((lhs-rhs)(PDEtools[dchange](y(x)=1/z(x),deq,[z(x)])));

zde := -diff(z(x),x)-x^2*z(x)^2+1

The numerical solution for z(x):

> zsol:= dsolve({zde, z(0)=1}, z(x), numeric, range = -2 .. 2);

Comparing the plots:

> with(plots):
   display([ odeplot(zsol, [x,1/z(x)], x=-2..2, colour=green),
        plot(solex, x=-2..2, discont=true)], view=[-2..2, -10 .. 10]);

It's unable to execute

seq(points(i), i = 1 .. t)

because t has not been assigned a value.

I think you want

> PLOT(ANIMATE(seq([POINTS(seq(points(i),i=1..n))],n=0..90)));

For making temporary assumptions, rather than permanent ones, it's usually better to use assuming rather than assume.

> limit(hommod, t=infinity) assuming d>0,Delta<=0;

                     0

> limit(hommod, t=infinity) assuming d>0, Delta > 0, 
    sqrt(Delta) < d;

                     0

> limit(hommod, t=infinity) assuming d>0, sqrt(Delta) > d, _C1>0;

                   infinity

Try curvarrows in the surfarro package of my Maple Advisor Database, www.math.ubc.ca/~israel/advisor

For example (after installing the database)

> with(surfarro):
   curve := [ cos(t), sin(2*t)];
   F:= diff(curve,t);
   T:= F/sqrt(F[1]^2+F[2]^2);
   curvarrows(T,curve,t=0..Pi, arrowcolour=blue,      
    scalefactor=1/2, arrowthickness=1,scaling=constrained,
    arrownum=4); 

[] and () have completely different uses.  If sol is a list, sol[2] is the second item of the list.  If that item happens to be a function of one variable, then sol[2](t) will evaluate it at some value t.