Your function definition doesn't work because the syntax of a function definition is

f:= (x,y)-> ...

not

f(x,y):= ....

But, instead of a function, it would be much easier to make f an Array, like this:

LL:= [[-11, -6, -1, 5, 5, 5, 5, 5], [-6, 1, 3, 6, 9, 9, 9, 9], [0, 5, 10, 13, 8, 6, 6, 6],
[9, 11, 12, 10, 10, 9, 9, 9], [9, 10, 11, 12, 12, 5, 5, 5], [11, 11, 11, 11, 10, 7, 7, 7],
[12, 12, 11, 13, 12, 7, 7, 7], [12, 12, 11, 11, 13, 8, 8, 8]]:
f:= 150+Array(0..7, 0..7, ListTools:-Transpose(LL)):
add(f[0,y], y= 0..7);
                              1236

Note that to save typing I made each original entry the difference of the number from 150. Then I added 150 to everything at the end.

Here's how (one way) to enter this system of PDEs into Maple. Like Preben, I need some clarification on the initial conditions.

restart:
alias(X= x(t,tau), Y= y(t,tau), Z= z(t,tau), Q= q(t,tau)):
xd,yd,zd,qd:= map(diff, [X,Y,Z,Q], t)[]:
xt,yt,zt:= map(diff, [X,Y,Z], tau)[]:
M:= (xt*(-Y-Z)+yt*(X+a*Y)+zt*(b-Z*(X-mu)))/(xt^2+yt^2+zt^2):
XM:= -Y-Z-xt*M:  YM:= X+a*Y-yt*M:  ZM:= b+Z*(X-mu)-zt*M:
pde1:= xd = XM+epsilon[x]*Q:
pde2:= yd = YM+epsilon[y]*Q:
pde3:= zd = ZM+epsilon[z]*Q:
pde4:= qd = a__x*XM+a__y*YM+a__z*ZM+beta*Q:

#Boundary conditions
BCs:= x(t,0)=x(t,2*Pi), y(t,0)=y(t,2*Pi),
     z(t,0)=z(t,2*Pi), q(t,0)=q(t,2*Pi),
     D[2](x)(t,0)=D[2](x)(t,2*Pi), D[2](y)(t,0)=D[2](y)(t,2*Pi),
     D[2](z)(t,0)=D[2](z)(t,2*Pi), D[2](q)(t,0)=D[2](q)(t,2*Pi)
:
#Initial conditions
# Need clarification on these.
ICs:= NULL:

#Parameters:
mu:= 1.85:  beta:= -2.23:  a:= 1/2:  b:= 3/4:
epsilon[x]:= -1/2:  epsilon[y]:= -1:  epsilon[z]:= -1/2:
a__x:= -1:  a__y:= 2:  a__z:= -1:

#Get the solution.
#The command is pdsolve, not dsolve, since these are PDEs, not ODEs.
Sol:= pdsolve({pde||(1..4)}, {BCs, ICs}, numeric):

Error, (in pdsolve/numeric) initial/boundary conditions must be defined at one or two points for each independent variable

How about

ex:= f(h*a[3]+x, y+h(a[3]-b[32])*k[1]+h*b[32]*k[2]):
mtaylor(ex, [x=0, y=0]);

Why not this?

xlist:=[seq(X1(i), i=1..10)]:
flist:=[seq(F(xlist[i],i), i=1..10)]:

The inner integrand has a singularity of order O(1/x^3) at x=0, so this can't be integrated over. This can be determined by looking at the leading term of the Laurent series.

f:= HankelH1(1, x*exp(I*Pi/4))^2/x;

series(f, x=0, 0);

Okay, here is a point plot in the complex plane of the roots of the numerator and denominator.

f:= 6*(20*z^5+137*z^4+450*z^3+850*z^2+900*z+420)
     /(300*z^6-1490*z^5+4197*z^4-7800*z^3+9600*z^2-7200*z+2520):
R1:= [fsolve(numer(f), complex)]:  R2:= [fsolve(denom(f), complex)]:
plot((L-> [Re,Im]~(L)) ~ ([R1,R2]), style= point, symbolsize= 16,
      legend= [numer,denom], labels= [Re,Im]);

The (only) command that you need to use is PDEtools:-dchange.

restart:
macro(E=E(x,y,z,t)):
PDE:= diff(E,x$2)+diff(E,y$2)+diff(E,z$2) = mu0*eps0*diff(E,t$2):
PDEtools:-dchange({x= r*cos(theta), z= r*sin(theta)}, PDE, simplify);

The error message indicates that broSet is not initialized. There needs to be a line somewhere such as

broSet:= {};

If that is not enough information for you to be able to solve the problem, then please post the whole code.

restart;
ode:= diff(y(t),t) = t*(1-0.3*t)-t*y(t)/(1+0.6*t):
ic:= y(0)=1:
Sol1:= dsolve({ode,ic}, numeric):
Sol2:= dsolve({ode,ic}, numeric, method= classical[foreuler], stepsize= 0.1):
P1:= plots:-odeplot(Sol1, [t,y(t)], t= 0..1, color= red, legend= rkf45):
P2:= plots:-odeplot(Sol2, [t,y(t)], t= 0..1, color= green, legend= euler):
plots:-display(P1,P2);

The problem is that your procedure WindowN does not know what to do when passed a symbolic value of Ci. Remember that arguments are evaluated before they are passed, so it's trying to evaluate WindowN(RP0,Ci,CangleWidth+Ci).

There are two solutions. The first is to use unevaluations quotes on the call to WindowN:

'WindowN(RP0,Ci,CangleWidth+Ci)'

The second, the one I'd prefer, is to have WindonN return unevaluated when appropriate. Include this line near the beginning  of WindowN:

if not args[2]::realcons then return 'procname(args)' end if;

You can (and should) substitute the actual second parameter name where I have args[2].

First, you misspelled restart.

The type command will not apply any simplifications before checking the type. So, for example,

x:= (sqrt(3)-sqrt(2))*(sqrt(3)+sqrt(2)):
type(x, integer);                         false

In other words, type only checks the properties that an expression has manifestly and explicitly. So, you need to simplify. Butthe simplify command is not good enough for this task. I did the following to your code:

p:=simplify((a+b+c)/2);
S:=expand(combine(sqrt(p*(p-a)*(p-b)*(p-c)), radical));
Now, setting N:= 4, I got 48 solutions in 5-8 minutes.

In Maple, an expression with an arrow (->) is equivalent to a procedure definition. So

f:= x-> ...f(x)...

defines f as a recursive procedure. All that is fine and legal in Maple.  Where you run into problems is when there is no base case to the recursion. Note that the error message was "too many levels of recursion", not "recursive assignment" which you would get if you attempted x:= x. The solution to your problem is to use unapply instead of the arrow, as Preben says.

I think that you have unfairly maligned Maple. Maple is no different in this respect from other languages that allow directly recursive procedures.

It seems that powsolve (and hence powseries) fails on this equation.

This dsolve solution doesn't explicitly derive one solution from the other, but it may be useful to you:

dsolve(myode, y(x), type= series);

Let V be the Vector of strings. Then do

sprintf(
     cat(
          "http://quandl.com/api/v1/multisets.csv?columns=",
          "%s," $ upperbound(V)-1, "%s"
     ),
     convert(V, list)[]
);

You need intial conditions for your system. I made up some completely arbitrarily for the code below. These are nonlinear ODEs, by the way---not that that distinction makes any difference in the code below.

sys:=
      diff(x(t), t) = 10*(y(t)-x(t)),
      diff(y(t), t) = x(t)*(28-x(t)) - y(t),
      diff(z(t), t) = x(t)*y(t) - 8*z(t)/3
;
ICs:= x(0)=1, y(0)=2, z(0)=3: #Initial conditions.
Sol:= dsolve({sys, ICs}, numeric):
plots:-odeplot(Sol, [x(t), y(t), z(t)], t= 0..100, numpoints= 10000);