I don't understand your line

  new defined parameters,  K= b/a = e/c = X / M =  Y/ R

X/M and Y/R are not likely to be constants.  In fact, if X(t)=  Xo- M(t),

X(t)/M(t) = Xo/M(t) - 1, which can't be a constant unless either Xo = 0 or M(t) is constant.

I don't think these differential equations will have closed-form solutions in general.  I think estimating values of the parameters from numerical data may be a difficult problem.

Hint: the sum of independent normal random variables is normal; the mean of the sum is the sum of the means, and the variance is the sum of the variances, of the individual random variables.

LSSolve in the Optimization package will do constrained least squares.  For example:

> Optimization[LSSolve]([seq(R[i], i=1..14)], sigma[1]=0..500, sigma[2] = 0..500, 
    sigma[3] = 0 .. 500, sigma[4] = 0 .. 500);

> plot3d([r,theta,f(r,theta)], r=0..1, theta=0..2*Pi, coords=cylindrical);

Or if you meant a two-dimensional plot of the curve f(r,theta) = 0,

> plots[implicitplot](f(r,theta), r = 0 .. 1, theta = 0 .. 2*Pi, coords=polar);

To my mind this is the most natural way to do it.

> map(`+`, A, a);

In Maple 13 you can also do it by

> A +~ a;

Actually there is a "Tutor" for this.  From the menu bar, choose:

Tools, Tutors, Linear Algebra, Linear System Solving..., and either Gaussian Elimination or Gauss Jordan Elimination.

But if you just want to solve the system, the simplest way is with the command

> solve({ 2*x-3*y+4*z =-12, x-2*y+z=-5, 3*x+y+2*z=1 });

From a theoretical point of view, the way to prove that the expected return is positive would be using some version of the Optimal Stopping Theorem for martingales and submartingales. 

Note that the exponential function is exp.  e is just another variable in Maple.

> cvalues:= [100, 200, 300, 400];
  plot([seq(900 + c*exp(t/2), c = cvalues)],t=0..5,
       labels = [t,p], legend=map2(`=`,c,cvalues));

I think you need to give fsolve a bit of help by cutting down the intervals where you're telling it to search.  Note that according to eq1, p <= pa + R*(d+pa*tan(beta)), and q/((1+alpha-alpha/cos(beta))) <= d+pa*tan(beta).

So try using

fsolve({eq1,eq2}, {p=0..pa + R*(d+pa*tan(beta)), 
  q=0..(1+alpha-alpha/cos(beta)) *(d+pa*tan(beta))})

Let's say P is a plot structure, which may involve the wrong colour.  For example:

> P:= plot(x^2, x=0..1):

To see what COLOUR (or COLOR) structure is there:

> indets(P, specfunc(anything, {COLOR, COLOUR}));

{COLOUR(RGB, 1.00000000, 0., 0.)}

Then you can change this colour to something else using subs.  For example, let's make it Navy:

> P:= subs(op(%) = COLOUR(RGB,op(ColorTools:-NameToRGB("Navy")/255.0)), P);

If the plot structure doesn't already have a COLOUR, things may be more complicated. 

Normally in PDE's, all the dependent variables are functions of the same independent variables.  So although

> pdsolve(sys);

does not work, you do get a result if you changed c(s) to c(s,n) and added a new equation diff(c(s,n),n) = 0.

However, this symbolic form of pdsolve does not allow boundary conditions, and does not produce a complete solution.  There is a numeric form of pdsolve, but I think it would require more in the way of initial and boundary conditions than what you have. 

Perhaps something like this?

> with(plots):
   P1:= plot3d(sin(sqrt(x^2+y^2)), x=0..4,y=0..4):
   P2:= plot3d([0,y,z], y=0..4, z=-1..1, colour=cyan, grid=[5,5]):
   P3:= plot3d([x,0,z], x=0..4, z=-1..1, colour=yellow, grid=[5,5]):
   P4:= plot3d([x,y,-1], x=0..4, y=0..4, colour=magenta, grid=[5,5]):
  display(P1,P2,P3,P4, axes=frame, orientation=[20,60]);

Do you really mean you want to solve for just the one variable y1s?  That's what your command means.  And that would only work if there was a value of that variable that made all 5 equations true (as functions of the other variables).  Since there is no such value, Maple quite correctly tells you there is no solution.  You could pick one equation and solve it for y1s, or you could solve your system of 5 equations for 5 variables.

Is this what you mean?

> diff(G(x(t,u), y(t,u)), t$2, u$3);

(the result is more than a bit complicated)

To enforce the constraint  0 <= y <= x, you could put in a new variable t with the constraint y = t*x and the interval
t = 0 .. 1.  Thus:

> fsolve({sin(x) + sin(y) = 1/2, cos(x) - cos(y) = 1, y = t*x}, 
       {x = 0 .. 2*Pi, y = 0 .. 2*Pi, t=0..1});

{t = .2947267864, x = 5.769236366, y = 1.700348494}