Given numeric values for B, G, then it *seems* as if one can obtain a numeric solution. Unfortunately different answers are obtained depending whether I use cartesians or polars :-(

For what it is worth (and treat the following with care, because at leasst one of the answers generated by the following must be wrong

   restart:
#
# Define integral
#
   P := (x, y) -> (1/2)*exp(-(1/2)*(x^2+G*y^2-2*B*x*y)/(-B^2+G))/(Pi*sqrt(-B^2+G));
#
# Change to polars
#
   P2:=(r, theta)->changecoords( P(x,y), [x,y], polar, [r,theta]);
#
# Try for symbolic solution in polars - takes ages
# doesn't seem to get anywhere - so comment out
#
#  int(P2(r,theta), theta=0..2*Pi, r=0..infinity);
#
# Try for numeric solution in polars
#
   Int(P2(r,theta), theta=0..2*Pi, r=0..infinity);
   evalf(eval(%,[G=2, B=-1]));
#
# Try the "same" integral in cartesians - OOOPs
# different answer
#
   Int(P(x,y), x=-infinity..infinity, y=-infinity..infinity);
   evalf(eval(%,[G=2, B=-1]));

Download tryStuff.mw

(slightly different) way

restart;
A:=Vector[column]( 4,
                                   [ J_K = Record(mu = 724.901557888305, sigma = 96.7437910529146),
                                     I_W = Record(mu = 775.098442111694, sigma = 96.7437910529198),
                                     K_J = Record(mu = 785.098442111694, sigma = 96.7437910529198),
                                     D_B = Record(mu = 764.901557888305, sigma = 96.7437910529146)
                                  ]
                               );
B:=sort(A, (i,j)->rhs(i):-mu > rhs(j):-mu );

Or you can sort 'in place' by using sort[inplace](...)

The error message you receive is because the output of your Maple command is very, very, VERY long. We are talking tens if not nundreds of screenfuls. You, generally, do not want to read this - except in very exceptional circumstances!!

Things you should consaider

1. Did you expect to get a really long output? If you did not, then something has gone wrong with the command which you executed. So fix the command!
2. If you did expect to get a really long output (but maybe it will get 'shorter' as you perform subsequent calculations), then suppress the long output by terminating the offending command with a colon character (ie ':') rather than a semicolon (ie ';') character
3. As Ronan has suggested, you can increase the limit on the allowed output - but you might want to think carefully about this. Depending on your specific installation and things like screen size/resolution etc, you might get about 80 characters per line, and 40 lines per screen - so that's about 3200 characters per screenful. (I could easily be off by a factor of two either way!) So 100000 characters is about 30 screenfuls - do you want/need 30 screenfuls? Increasing the output limit (as Ronan suggested) will give you more 'screenfuls'. Do you want that?????

because they are not independent. You can determine this pretty much by 'eyeball' - swapping x2 and x3 in the defined expressions.

Alternatively you can use

  restart;
#
# Define equations
#
   eq1:= 1752*x1+384*x2+384*x3+ 72*x4=0:
   eq2:=  384*x1+366*x2+ 42*x3-144*x4=0:
   eq3:=  384*x1+ 42*x2+366*x3-144*x4=0:
   eq4:=   72*x1-144*x2-144*x3+216*x4=0:
#
# Solve equations
#
   sol:=solve([eq1,eq2,eq3,eq4],[x1,x2,x3,x4]);
#
# The above 'solution' states that x3=x3 - in
# other words x3 is absolutely arbitrary. The
# system will be solved if one sets x2=x3
# (with x3 arbitrary), x1=-x3/2 (with x3
# arbitrary) and x4=3*x3/2 (with x3 arbitrary)
#
# As a check substitute suggested values for x1
# (ie sol[1][1]),  x2 (ie sol[1][2]) and x4
# (ie sol[1][4]) in [eq1, eq2, eq3, eq4]
#
#
   s2:=subs( [sol[1][1], sol[1][2], sol[1][4]],
                     [eq1, eq2, eq3, eq4]
          );

Have to admit that I like rlopez' solutiion.

But just to illustrate that there is generally more than one way to do this - here are another two

restart;PP := .8707945038*exp(-50.00000000*(m-0.842e-1)^2+(2.745342070*(m-0.842e-1))*(a-2.3722)-.1046792095*(a-2.3722)^2):
plt1:=plot3d(PP, a=-0.5..0.5, m=-0.5..0.5,color=red, style=surface, transparency=0.5):
plt2:=plot3d( [x, -0.2, eval(PP,[a=x, m=-.2])], x=-0.5..0.5, y=-0.5..0.5, thickness=5):
plots:-display( [plt1, plt2], axes=boxed, view=[-0.5..0.5, -0.5..0.5, -0.5.0..1.0] );

or

restart;
with(geom3d):
PP := .8707945038*exp(-50.00000000*(m-0.842e-1)^2+(2.745342070*(m-0.842e-1))*(a-2.3722)-.1046792095*(a-2.3722)^2):
plt1:=plot3d(PP, a=-0.5..0.5, m=-0.5..0.5,color=red, style=surface, transparency=0.5):
plt2:=plot3d( [x, -0.2, eval(PP,[a=x, m=-.2])], x=-5..5, y=-5..5, thickness=5):
point(p1, [-5, -0.2, eval(PP, [a=-5, m=-0.2])]):
point(p2, [ 0, -0.2, eval(PP, [a= 0, m=-0.2])]):
point(p3, [ 5, -0.2, eval(PP, [a= 5, m=-0.2])]):
plane(pl1, [p1,p2,p3]):
plt3:=draw(pl1, style=surface, color=cyan, transparency=0.5):
plots:-display( [plt1, plt2, plt3], axes=boxed, view=[-0.5..0.5, -0.5..0.5, 0.0..1] );

I just want to you to do some *serious* thinking about your question.

1. It is trivial to solve your ODE system for any given value of 'w'
2. It is trivial to plot all the dependent functions ( ie x(t), y(t), z(t) ) in your ODE system, over any range of the independent variable 't' (and do this for any value of the parameter 'w')
3. It is trivial to extract the value of all dependent variables ( ie x(t), y(t), z(t) ) in your ODE system, at any given vaue of the independent variable 't' (and do this for any value of the parameter 'w')
4. The problem I have is that nothing changes very much! x(t), y(t), z(t) are all pretty much constant in the variable 't'. And they don't vary very much as the parameter 'w' is varied

The following code will alow you to experiment with the above just change the assignment w:=100 to any numeric value and examine the solutions. Now do you really want to spend a lot of time, genrating a lot of solutions which don't vary significantly???

  restart;
#
# Define Parameters. NB parameter w is set to 100
# although setting it to 1000 or even 10000 doesn't
# make much difference
#
   alpha:=.5:
   gamma0:=8.82*1e10:
   mu0:=4*Pi*1e-7:  
   gammA:=gamma0/(1+alpha^2):
   alphaPrime:=mu0*gammA*alpha:
   gamma1:=mu0*gammA:
   m0:=mu0*1e5:
   wl:=0.08:
   w:=100:
#
# Define ODEs and initial conditions
#
   dsys:={ diff(x(t), t) =  alphaPrime*(z(t)^2*wl*cos(w*t)/gamma1+y(t)^2*wl*cos(w*t)/gamma1),
                 diff(y(t), t) = -z(t)*wl*cos(w*t)-alphaPrime*y(t)*wl*cos(w*t)*x(t)/gamma1,
                 diff(z(t), t) =  y(t)*wl*cos(w*t)-alphaPrime*z(t)*wl*cos(w*t)*x(t)/gamma1,
                 x(0)=.01,
                 y(0)=1,
                 z(0)=0.1
             }:
#
# Solve ODE system
#
   dsn1 := dsolve(dsys, numeric, maxfun=0):
#
# Get values of all solution variables at t=1
# and t=100, just as a check. Note that these
# values do not change very much
#
   dsn1(1);
   dsn1(100);
#
# Plot values of all dependent functions for the range
# t=0..10
#
   plots:-odeplot( dsn1,
                              [ [t, x(t)],
                                [t, y(t)],
                               [t, z(t)]
                             ],
                             t=0..100,
                             color=[red, blue, green]
                          );

Works for me!

I suppose it *might* have been some temporary issue with Maple servers?

For example, your issue with fsolve() can be achieved by using the setmaple() command

   setmaple('Digits', 50)
   syms x
   maple( 'fsolve', 'x^3-4*x^2+3*x+6')
   ans = -0.84546609143593277372176028502548481300690617199393

It is necessary to use setmaple() in order to change settings in the current Maple session, which will be available for subsequent commands. For example

   setmaple('a', 'cos(x)')
   getmaple('a')
   maple( '2*a')

assigns the variable 'a' to cos(x) within the Maple session, so the subsequent two commands will return

   ans = cos(x)
   ans = 2 cos(x)

as you would expect.

One 'feature' (which is a trap for the unwary) is that some Maple commands are available without the maple() wrapper, and some are not. So

syms x

int(2*cos(x),x)

diff(2*cos(x), x)

will all work. However

dsolve( diff(f(x),x)=cos(x) )

will fail, because one has to use

maple( 'dsolve', 'diff(f(x),x)=cos(x)')

The list of Maple commands which will work in Matlab without the the maple() wrapper can be examined, by accessing the Maple Toolbox help within Matlab. If you fire up the Matlab help, then the left hand pane *should* contain two entries, "My Products" and "Supplemental Software". Under the latter should be "Maple Toolbox" - just click it and you now have access to all the relevant help, including a few examples. This 'help' is of a similar standard to the Maple help within Maple, or indeed the Matlab help within Matlab

You prettyprint comment is a bit ambiguous. I think you have to realize that although you are using the  Maple engine for calculation, you are using the Matlab interface for entering input and displaying output - so you have to live with the limits of this Matlab interface. Roughly speaking this corresponts to the interface(prettyprint=1) setting withiin Maple

f:=(x,y)->(x^3+y^3)^(1/3);
g:=(x,y)->diff( f(x,y),x);
g(x,y);
eval( g(x,y), x=1);

You need to resolve the following, before your question makes any sense at all

1. You seem to be using E[s], E[c], G[c] and  E__c, E__s, G__c, more or less interchangeably. Don't. Pick one and stick to it - I would recommend the latter for your problem

2. You use the parameter 'l' (that's a lower-case L) although this parameter is never defined. L[0], L[1], L[2], L, have all been defined. You just have to determine what you mean when you use 'l'

3. You use the parameter Zeta (note the case) in a few places, although this has never been defined. It might be  zeta__c or zeta__s??

4. Your final system contains the dependent functions f3(x), u__ob(x, t), u__ot(x, t), w__b(x, t), w__t(x, t). Note that four of these are functions of two variables - so you have a system of partial differential equations. This will never be solved using dsolve(), although pdsolve() might work if all other syntax errors were fixed

5.  

   I can't guarantee to have located all syntax errors

    

    

    

    

because I can't subsequently incorporate the bounday condition. Bu the following will provide a 'solution' for your PDE

  restart;
#
# Define PDE
#
  PDE:=diff(c(x,t),t)+c(x,t)*diff(c(x,t),x)=0;
#
# Look for a separable solution of the form
# c(x,t)=f(x)*g(t), and 'build' the result
#
  sol1:=simplify(pdsolve(PDE, HINT=`*`, build));
#
# Verify this solution by substituting in the
# the original PDE - should get 0=0
#
  subs[eval](sol1, PDE);

But the following seems to provide the same answer(s) and is a lot shorter

coeffs.mw

You seem to want a numeric rather than a symbolic solution so use fsolve(). It runs more or less instantaneoudly, as in

Download numerSols.mw

As in

ygraph1 := -.736312023696564122*exp(2.26140104440167664*10^5*tt)-.591826613918776445*exp(28994.5376895644186*tt)+.328002839648234568*o*exp(13767.7178702679158*tt);
ygraph2 := -.591859486202007235*exp(2.26140104440167664*10^5*tt)+.328381376616263988*exp(28994.5376895644186*tt)-.736116852194203974*o*exp(13767.7178702679158*tt);
ygraph3 := -.327943520064913564*exp(2.26140104440167664*10^5*tt)+.736143281263262450*exp(28994.5376895644186*tt)+.592069351595225779*o*exp(13767.7178702679158*tt);
 

You may have to scroll right to see the offending character(s)