Do you want to plot it, or convert it to rectangular?  Do you realize that  your r will be negative, or do you really mean r = sin(theta-Pi)?

Plotting is the simplest:

> plot(sin(theta)-Pi, theta = 0 .. 2*Pi, coords=polar);

To convert an expression in polar coordinates (r, theta) to rectangular coordinates:

> eval(expr, {r = sqrt(x^2+y^2), theta = arctan(y,x)});

For dsolve, a system should be a set or list, not just a sequence of equations.  Thus

> dsolve({sysode});

It looks to me like the problem is your use of exp as a local variable.
That wouldn't be a problem, except that you also want to use the exponential function, so your code contains

exp:=-diff(Asy(t),t)+(1-exp(-0.8*t))+.7;

which makes exp impossible to evaluate.

If you really want to call your local variable exp, you could use :-exp for the exponential function, but why not just use a different name for the local variable?

> with(plots):
  complexplot((-1/2)^q+(-1/2)^(1-q), q =  0.25 .. 0.5);

Actually complexplot is nothing more than a wrapper for Alex's solution.

First you should state the assumptions (unfortunately, these are often not explicitly stated in this sort of problem).  A person's birthday is, supposedly, equally likely to be on any day of the year (for convenience, February 29 is usually ignored), and different people's birthdays are independent. 

Next, here's a hint: it's often easier to consider the complement of an event rather than the event itself.

See www.mapleprimes.com/forum/kronecker/tensor-products and www.mapleprimes.com/blog/jacquesc/maple-12-wish-list

It's impossible to map a region of the plane onto a region of a sphere isometrically (i.e. without changing distances, as measured along the surface).  This is because the Gaussian curvature of the sphere is positive, while that of the plane is 0.  See e.g.
en.wikipedia.org/wiki/Curvature#2_dimensions:_Curvature_of_surfaces
 

This is an interesting example. 

It looks like Du should be (ke+ce*s) times the 3 x 3 identity matrix.

> Probe:= map(t -> (t = rand()),  indets(Pt) minus {ke, ce, s}));

 >  (subs(Probe,G21) . (subs(Probe, G11) - subs(Probe,Pt))^(-1) . subs(Probe,G12) - subs(Probe,G22))^(-1);

  [  ke + s ce     0         0      ]
  [     0       ke + s ce    0      ]
  [     0          0      ke + s ce ]

Presumably one could prove this by a sufficient number of different substitutions,
using bounds on the degrees of the numerators and denominators of the rational
functions, but a few repetitions is pretty convincing evidence.

First you need to get those digits of Pi in a string.

> PiDigits:= StringTools[Delete](convert(
        evalf(Pi,100000),string),2..2):

Now for a search procedure:

> SearchPi:= proc(n::nonnegint)
       SearchText(convert(n,string), PiDigits)
   end proc;

For example:

> SearchPi(1592653);

                     4

If the string is not found, it returns 0.

For your type:

> TypeTools[AddType](mytype, 
     t -> evalb(t::nonnegint and SearchPi(t) > 0));

As long as Digits > evalhf(Digits) (which is 14 on my system), the plotting routines should not use hardware floats.  Setting infolevel[plot] to 2 or more should confirm this.  For example:

> q2:= (x) -> x*exp(x^2)*(1-erf(x)):
   Digits := 50:
  infolevel[plot] := 2:  
  plot(q2(x), x = 0 .. 10);

plot/adaptive:   using evalf
plot/adaptive:   produced    55    output segments
plot/adaptive:   using    55    function evaluations

The entries of your matrices are huge expressions in many variables, and it's not too surprising that Maple has trouble computing them.  But with numbers substituted for most of the parameters, everything should be much easier.  My suggestion is that you avoid as much as possible actually calculating things such as Du until putting in those numbers.  For example, rather than having an actual matrix Du, have a procedure that substitutes values for the parameters and then calculates Du.

The name "singularperturbation" should give a clue that there are likely to be severe problems here with any numerical method.  That seems to be happening here.
I tried this:

> ivp := { 0.01 * diff(y(t),t$2) - y(t)*diff(y(t),t) = 0, y(1) = -1, D(y)(1) = m } ;
   dsolve(ivp);

y(t) = 1/10*tan(1/2*(10*t*(m-50)^(1/2)-2^(1/2)*arctan(5*2^(1/2)/(m-50)^(1/2))-10*(m-50)^(1/2))*2^(1/2))*(m-50)^(1/2)*2^(1/2)

> ys:= rhs(%);
   ym := eval(ys, t=-1);
   plot(ym, m=-1 .. 1);

It seems that to a very good approximation, m=0 would make y(-1) = 1.
That can't be exactly right, because it's clear that the solution of the
differential equation with y(1) = -1, y'(1) = 0 is the constant y = -1.

What have you tried so far, and where are you stuck?  You should not be asking us to do your homework for you, you should be asking for hints on specific points that you have difficulty with.

I don't know if it's worth the trouble, but you can try this:

> eq:= 7*x = 14;
   isolate(subs(map(t -> (t = ``(t)),indets(eq,integer)),eq),x);

x = ``(14)/``(7)

> expand(%);

x=2
  

I presume first-order conditions are the Karush-Kuhn-Tucker conditions, which are necessary (but not sufficient) conditions for a maximum or minimum [ assuming the constraint qualifications ] involving the first partial derivatives of the objective and constraints.  For example, for an interior point of the feasible region the first-order conditions say the gradient is 0.