There were quite a few syntax issues. 

ode.mw

The best is to use piecewise, see for details maple help

restart:with(plots):
p1 := proc (x, y) piecewise(x^2+y^2 <= 1 ,x*y-y^2,0) end proc;
plot3d(p1(x,y), x = -1 .. 1, y = -1 .. 1);

To avoid the error you are getting, just use single quotes

p1 := proc (x, y) if x^2+y^2 <= 1 then x*y-y^2 else 0 end if end proc;

plot3d('p1(x,y)', x = -1 .. 1, y = -1 .. 1);
 

As mentioned by Carl Love, your question is unclear. But if you want to see the ouput of the two loops then you need to use printlevel:=3 before the loops. 

The so-called analytical solution you proposed is not satisfying the ode,

Eq:=diff(w(x), x$6)+diff(w(x), x$4)+diff(w(x), x$2)-w(x) = -90;

sol1:=w(x)=C1*sinh(x)+C2*cosh(x)+C3*sin(x)+C4*cos(x)+C5*sin(x)+C6*cos(x)+90;

odetest(sol1, Eq);
 

To verify a solution to an ode, you can use odetest to do so.

Download 1.mw

As for as I know, the dsolve has already given you the general solution in compact form.

Download 1.mw

You can also opt for the continuation option.

Download cErrore.mw

In Maple 2016, I have executed your sheet without any problem.

Question.mw 

You can use the many options available  with the 'contourplot' to improve the quaility.

Here I used the `contours=s` and `grid`

contourplot(-(1/2)*y^2-(1/2)*x^2-(1-.3)/sqrt((x+.3)^2+y^2)+((-1)*.3)/sqrt((x-1+.3)^2+y^2), x = -1.5 .. 1.5, y = -1.5 .. 1.5,grid=[60,60],contours=120);

You can get the 2d look from the 3d just by using orientation option

contourplot3d(-(1/2)*y^2-(1/2)*x^2-(1-.3)/sqrt((x+.3)^2+y^2)+((-1)*.3)/sqrt((x-1+.3)^2+y^2), x = -1.5 .. 1.5, y = -1.5 .. 1.5, view = -2 .. -1.3, axes = boxed,style=patchnogrid, grid=[160,160], orientation=[90,360], lightmodel=light4, shading=zhue,  transparency=0.3);
 

As mentioned by @ Kitonum, there is no built-in Maple command to find Fourier series of a function.

But you can tell Maple to help you avoid all the integration and summation.

To find a Fourier series of a function, first we need to find the Euler-Fourier formulas (a_m and b_n) and then sum it.  This can be done by writting three procedures like this

Download FS.mw

This not my own idea. I got it somewhere on the internet.

restart:with(plots):

x:=2*m*p/(m^2+n^2+p^2);

y:=2*n*p/(m^2+n^2+p^2);

z:=(p^2-(m^2+n^2))/(m^2+n^2+p^2);

p:=2:

plot3d([x,y,z],n=-1..1,m=1..1,grid=[30,30],scaling=constrained,axes=boxed);

f:=x^2+y^2=z^2;

implicitplot(f,n=-1..1,m=-1..1);

Just type it down in maple

restart:
B[d]:=(-d*w[1]+w[3])/(((-2*w[1]+w[3])^2)-((-d*w[1]+w[3])^2)):

B[o]:=(-2*w[1]+w[3])/(((-2*w[1]+w[3])^2)-((-d*w[1]+w[3])^2)):

A[o]:=w[1]*(alpha[o]-c[o]-t[o])+2*w[2]*e[o]:

A[d]:=w[1]*(alpha[d]-c[d]-t[d])+2*w[2]*e[d]:

q[o]:=B[o]*A[o]-B[d]*A[d]:

q[d]:=B[o]*A[d]-B[d]*A[o]:

U[d]:=w[1]*(q[d]*(alpha[d]-q[d]-d*q[o]-c[d]-t[d])-C)+w[2]*(e[d]*q[d]+e[o]*q[o])+w[3]*((1/2)*(q[d]+q[o])^2);

restart:

for i to 2 by 1 do
n[i] := (1/2)*i
end do;

for j to 2 by 1 do
eq[j] := diff(f(eta), eta, eta)+(1/(4*eta^2)+(3.*n[j]+1)/(n[j]+1)*(NU/(4*a*eta)-NU*(eta/a)^((n[j]+1)/(2.0*n[j]))/(4*a*eta)))*f(eta) = 0:
end do;

for j to 2 do
dsolve(eq[j]);
end do;

As @Rouben mentioned, you need to post the ODEs you want an answer to. 

Checking the maple help page for dsolve, I found an example which might help

 Define a system of ODEs.

sys_ode := diff(y(t),t) = x(t), diff(x(t),t) = -x(t);

If the unknowns are not specified, all differentiated indeterminate functions in the system are treated as the unknowns of the problem.

dsolve([sys_ode]);

{x(t) = _C2 exp(-t), y(t) = -_C2 exp(-t) + _C1}

Define initial conditions.

ics := x(0)=1, y(1)=0;
ics := x(0) = 1, y(1) = 0

Solve the system of ODEs subject to the initial conditions ics.

dsolve([sys_ode, ics]);

You can then extend for a system of second order ODE's fairly easy. While trying for that if maple is unable to solve analytically your unseen problem then you can use numeric.

GoodLuck

Download FDM_Dbvp_sheet_1.mw