Your system doesn't have any basin of attraction: the equilibrium points are all centres and saddles. 
The trajectories lie on curves y^2 - 2*k*cos(x) = constant.  Here's a typical look at the phase portrait (taking k = -1):
 

> DEtools[DEplot]([D(x)(t)=y(t),D(y)(t)=sin(x(t))],[x(t),y(t)],
    t=-15..15,x=-6.5 ..6.5,y=-2.6 .. 2.6, 
    [[x(0)=0,y(0)=0.01], [x(0)=0,y(0)=-0.01],
     [x(0)=3,y(0)=2.5],  [x(0)=3,y(0)=-2.5],
     [x(0)=-3,y(0)=1],   [x(0)=3,y(0)=1]],
    linecolour=blue,stepsize=0.01);

The main problem is that you have two independent variables, x and t, which would make it a partial differential equation.
However, since there are no derivatives with respect to x, it's best to consider x as a parameter.  Now it looks to me like  you have a
constant-coefficient linear system.  You should be able to get a symbolic solution, not just a numerical one, using linear algebra.  For symbolic solutions, you should work with the complex system, not break it down into real and imaginary parts.  Something like this I think:

> sys := [diff(Bx(t), t) = I * x*Bx(t) + I * Ux(t),
          diff(By(t), t) = I * x * By(t) + I * Uy(t) - 3 * Bx(t)/(2*omega),
          diff(Bz(t), t) = I * x * Bz(t) + I * Uz(t),
          diff(Ux(t), t) = x * Ux(t)  + I * Bx(t) - 2 * I * Uy(t)/omega, 
          diff(Uy(t), t) = I * x * Ux(t) + I * Bx(t) - Ux(t)/(2*omega),
          diff(Uz(t), t) = I * x * Ux(t) + I * Bx(t)];

> with(LinearAlgebra):
   A, b := GenerateMatrix(map(rhs,sys), [Bx(t),By(t),Bz(t),Ux(t),Uy(t),Uz(t)]);

                  [x I ,    0 ,    0 ,   I ,   0 ,   0]
                  [                                   ]
                  [     3                             ]
                  [- ------- ,  x I ,  0 ,  0 ,  I , 0]  [0]
                  [  2 omega                          ]  [ ]
                  [                                   ]  [0]
                  [0 ,    0 ,    x I ,   0 ,   0 ,   I]  [ ]
                  [                                   ]  [0]
          A, b := [                        -2 I       ], [ ]
                  [I ,   0 ,   0 ,   x ,   ----- ,   0]  [0]
                  [                        omega      ]  [ ]
                  [                                   ]  [0]
                  [                 1                 ]  [ ]
                  [I , 0 , 0 , - ------- + x I , 0 , 0]  [0]
                  [              2 omega              ]
                  [                                   ]
                  [I ,    0 ,    0 ,   x I ,   0 ,   0]

> P:= factor(CharacteristicPolynomial(A,r));

P := -(x+r*I+2*I*omega-2*I*omega*x^2+2*x*r*omega-r*omega^2-I*omega^2*x^2*r+r^2*omega^2*x*I+r^2*omega^2*x-r^3*omega^2)*r*(-r+x*I)^2/omega^2

Unfortunately the cubic factor is going to cause trouble, because Maple will try to use the very complicated solution of the cubic in radicals.

There are two ways to concatenate strings, || and cat.  They differ in their evaluation rules, but that's not too important here.  Note also that name should not be used for a variable name, as it has an important meaning within Maple..
 

> myproc:= proc(n)
     local i;
     i:= 10;
     save i, cat("i_",n);
     ExportMatrix(cat("Matrix_",n), <<1,2>|<3,4>>);
     end proc;

I suspect you'll want to look at the real and imaginary parts.  Replace any variable z(t) that will be complex with x(t) + I*y(t),
and use evalc(Re(...)) and evalc(Im(...)) to take real and imaginary parts of your equations, giving you a system of DAE's for real variables.

What's your question?  I know what a simplex is, but I don't understand exactly what you want to do.

Comments should start with #, not *.  Other than that, it works.  I would use symbol=point in the pointplot, though, and put the unassign('r') in the Bifurcation procedure.  Also, you don't need to use display, just

> P1; P2;

complexplot is for plotting a parametric curve in the complex plane.  You don't have that.  You might try implicitplot (after splitting z into real and imaginary parts).  Thus:

> with(plots):
  eq:= abs(x+I*y-a) - abs(x+I*y+a)=2*c;
  implicitplot(eval(eq, {a=I,c=1/2}), x=-2..2, y=-2..2);

In principle you could program a Go-player in Maple, but in practice it would be much too slow (except maybe for prototyping).  For something this computation-intensive, you really want a low-level compiled language such as C,

Or, for something more like the Mathematica version:

> s:= FileTools[Text][ReadFile]("C:\\Documents and Settings\\Don Taylor\\Desktop\\sample.srt");
  with(StringTools):  
  for line in Split(s,"\n") do 
     if not RegMatch("[0-9][0-9]:[0-9][0-9]:[0-9][0-9].*", line) 
      then printf(cat(line, "\n")) 
     end if 
  end do;

Your fsolve call makes no sense as written.  The first argument is a set of two procedures with four arguments each.  And then you have  a second argument which contains two constants (actually one because the two are equal), an equation and an inequality.  Note that there is no connection between the formal parameters e1, e2, As1, As2 of the procedures and the global variables of the same names.

Exactly what equations are you trying to solve, for what variables?  Note that you need as many equations as variables.

When y = 0, your expression x*sqrt(1-y^2) + y*sqrt(1-x^2) has real values for all real x.  When y is not zero, it only has real values for -1 <= x <= 1.  This shows up here, as plot3d joins a point with x > 1, y = 0 to two neighbouring points with x=1, forming a little triangle that sticks out of the rest of the plot.

First of all, it's better to use Matrix rather than matrix.

> c:= Matrix([[a],[b]]);

> c^%T;

The system of DE's probably has no closed-form solution.  It can, of course, be solved numerically, given initial values and values of the parameters.  And then you can try LSSolve to fit parameter values to the data.  It may be tricky, though.  One problem is the rather nasty singularity at (X,Y)= (0,0), which might be reached in finite time along certain trajectories.

Doug is right: high-order polynomial interpolation is a numerical minefield.  Spline interpolation might work better.

But there's another way to solve your problem, by looking at a bigger DE system.

Consider, slightly more generally, the initial value problem

Y' = F(Y), Y(to) = Yo.

You want to find the derivative of the solution wrt a parameter p, where to and Yo are functions of p.  I'll write Y = Y(t, p).
Let Z(t,p) = (d/dp) Y(t,p).

If Y(to(p),p) = Yo(p), then Y(to(p+dp), p+dp) = Yo(p+dp) = Yo(p) + Yo'(p) dp + O(dp^2)

but according to the differential equation Y(to(p+dp), p+dp) = Y(to(p), p+dp) + F(Yo(p)) to'(p) dp + O(dp^2)
so Y(to(p), p+dp) = Yo(p) + (Yo'(p) - F(Y(to(p),p)) to'(p)) dp
i.e. Z(to(p), p) = Yo'(p) - F(Yo(p)) to'(p)
while (d/dt) Z(t,p) = (d/dp) (d/dt) Y(t,p) = (d/dp) F(Y(t,p)) = (grad F)(Y(t,p)) . Z(t,p).

So we have an enlarged initial value problem for both Y and Z, which can be solved numerically.

In your example, Y = [x, y, y'], F(x, y, v) = [ 1, v, (c + d*x^2)*y], to(p) = p, Yo(p) = [p, exp(-p), 1/2]:
You end up with this system (with yp = dy/dp and vp = dv/dp):

> des:= {D(y)(x)=v(x),D(v)(x)=(c+d*x^2)*y(x),
         D(yp)(x)=vp(x),D(vp)(x)=(c+d*x^2)*yp(x)};
  ics:= {{y(p)=exp(-p),v(p)=1/2, yp(p)=-1/2-exp(-p), vp(p)=-(c+d*p^2)*exp(-p)};

For example, with p = 1/6, c = 10, d = 1,

> Sol := dsolve(eval(des union ics, {p=1/6, c=10, d=1}), numeric);

At x = 2 we have

> Sol(2);

[x = 2., v(x) = 842.878276886481899, vp(x) = -3370.96972515993458, y(x) = 229.537728891821984, yp(x) = -917.995212516255606]

so that (d/dp) y(2) = -917.995212516255606.

Unfortunately densityplot ignores the option coords=polar.  However, there is a work-around.  You can try something like this.

> with(plots):
  f:= (r, theta) -> r^2*cos(2*theta);
  display(changecoords(densityplot(f(r,theta), r=0..1, theta = 0 .. 2*Pi, style=patchnogrid, 
     colorstyle=HUE), polar), axes = box, orientation = [270,0], labels = [x,y,``]);