For phi0=0 we can get an answer, although it takes a few seconds. It is integrating.
restart;
f:=phi->exp(-(phi-phi0)^2/2/sigma^2);
u:=eval(f(phi),phi0=0);
inttrans[fourier](u,phi,omega) assuming sigma>0;

From this you get the result for phi0<>0 by multiplying the above by exp(-I*omega*phi0) (since it just requires a simple change of variable).
Or do:
inttrans[fourier](g(phi-phi0),phi,omega) assuming phi0::real;

In the second case you have plot instead of plot3d, but your expression has two variables r and phi, thus you need plot3d.
The coordinate systems you can use are described in the help pages
? coords
Spherical coordinates are surely not what you want but the graph looks funny. Notice that I wrapped the expression in ln because of the huge numbers appearing:
plot3d(ln(abs(psi(r,phi))^2),r=-0..5, phi= 0..2*Pi, coords = spherical);

Addendum.

Maybe in the first case you want z_cylindrical coordinates:

addcoords(z_cylindrical,[z,r,theta],[z,r*cos(theta),r*sin(theta)]);
plot3d(r^2,r=0..5,phi=0..2*Pi,coords=z_cylindrical); #The order of the r and theta ranges matter!
Compare with
plot3d(r^2,r=0..5,phi=0..2*Pi,coords=cylindrical); #!!

You can use mtaylor instead. It is meant for functions of several variables, but works for a single variable too:
ex := mtaylor(exp(x), x = 0, 5);

As Carl is saying, the issue is really clarifying what it is you have as input, and what you want as output.
Assuming as I do that you want to input a function (procedure) and expect a function (procedure) as output, then you can do:

f:= proc(P) local x; unapply((x-1)*diff(P(x),x)-x*P(1),x) end proc;
#Notice here that it is assumed that P(x) makes sense, i.e. P is a function:
#Further, the use of unapply makes the output a function too.
f(q); #q not assigned
f(x->x);
f(1); #Maple accepts that a number can work as a constant function 45(x) = 45 for all x
f(x->1); # same result as above
f(sin);
f(sin)(t);
p := x-> 1+x+x^2;
f(p);
expand(f(p)(t));

sum evaluates its arguments before it does anything else (the normal way in Maple).
Thus f(4) first results in sum(x[newc],c=1..4), then 4*x[newc].

You can do
f:=C->sum('x[cat(new,c)]',c=1..C);
or use add:

f:=C->add(x[cat(new,c)],c=1..C);

To illustrate the point about evaluation:

f:=C->Sum(x[cat(new,c)],c=1..C); #Capital S
f(4);
value(%);

You can play around with this:

pp:=pointplot3d(Tab,color = magenta,axes = boxed,symbol = solidsphere,symbolsize = 50);
eps:=0.1:
nms:=[seq(seq(seq([Zx[i]-eps,Zy[j]-eps,Zt[k]-eps,u[i,j]^k],i=1..N+1),j=1..N+1),k=1..N+1)];
tp:=textplot3d(nms, 'align'='left', 'axes'='boxed');
display(tp,pp);

When you write a procedure p, you can at appropriate places insert lines like
userinfo(2,'procname',"Print some information");

If you outside the procedure do
infolevel[p]:=2; #or higher
then a call to p, say p(x), will result in the printing of "Print some information" besides p doing what it was made to do.
See
?userinfo
for a description and an example.
There you also find the statement:

Throughout the Maple library userinfo statements have been used with the following conventions:
  Level 1: reserved for information that the user must be told.
  Level 2,3: general information, including technique or algorithm being used.
Level 4,5: more detailed information about how the problem is being solved

With z=x+y your expression is

1/2*( -1 + sqrt(-3*z^2)/z)

Thus it is really a question of what the value is of w=sqrt(-z^2)/z.

Write z = r exp(I*theta) where r>0 and theta is real.
Then
-z^2=r^2*exp( I*(Pi+2*theta)) and
1/z= 1/r*exp(-I*theta).
Thus

w=sqrt(-z^2)/z = exp(-I*theta)*sqrt(exp(I*(Pi+2*theta)).
But Maple gives you the principal squareroot.
Suppose first that  -Pi<theta<=0, then -Pi<Pi+2*theta <= Pi,
so -Pi/2 < 1/2*(Pi+2*theta) <=Pi/2.
Thus in this case sqrt(exp(I*(Pi+2*theta)) = exp(1/2*I*(Pi+2*theta)) = exp(I*(Pi/2+theta))
and so w=exp(I*Pi/2 ) = I.

If on the other hand
0<theta<=Pi, then Pi<Pi+2*theta <= 3*Pi, so here the principal square root is
sqrt(exp(I*(Pi+2*theta)) = sqrt(exp(I*(Pi+2*theta-2*Pi)) = exp(I*(-Pi/2+theta)),
and so w=exp(-I*Pi/2 ) = -I.
###This was edited.

Two cases to illustrate this in Maple:
w:=sqrt(-z^2)/z;
evalc(eval(w,z=2-3*I));
evalc(eval(w,z=2+3*I));

Try
interface(warnlevel=0);

Try this:

N:=4.5: #Works in Maple 18
bcs1 := f(0) = 0, (D(f))(0) = 1, D(f)(N) = 0,F(N) = 0,G(N) = -f(N),H(N) = n,theta(0) = 1,theta(N) = 0,theta1(N) = 0;

print doesn't produce output, it only prints.
Try

f:=x->print(x^2);
a:=f(8);
a;
Just remove print:
f:=x->x^2;

Since I only have the image you presented to work with and not the code to copy and paste, I can only say that I definitely think you should change your dots (.) to multiplication signs (*) all over.

Multiply by a test function phi with compact support. Integrate both sides from -infinity to infinity.
The part having u''(x)*phi(x) you integrate by parts twice.
Then you will have something like
Int( (diff(phi(x),x,x)-(1+x^2)*phi(x) )*u(x),x=-infinity..infinity)=Int(phi(x),x=-infinity..infinity);

Now consider the ode:
diff(phi(x),x,x)-(1+x^2)*phi(x)=psi(x);
where psi is a test function with compact support in (-K,K). Find phi s.t phi(-K)=0 D(phi)(-K)=0.
This phi will also be a test function with compact support. Thus the inititial step holds for that phi.
That problem can be solved by Maple. The result for phi is given in terms of an integral of psi.
The shape is
phi(x)=Int(f(x,z)*psi(z),z=-K..x);
Now you have from the initial step that
Int(psi(x)*u(x),x=-K..K)=Int( Int(f(x,z)*psi(z),z=-K..x), x=-K..K);
Change the order of integration for the iterated integral on the right to get
Int(psi(z)* Int(f(x,z),x=z..K), z=-K..K);
Change the names for easy comparison with Int(psi(x)*u(x),x=-K..K).
Thus you get:
Int(psi(x)*(u(x)-Int(f(z,x),z=x..K),x=-K..K) = 0.

This is very sketchy and not finished.

Searching in MaplePrimes led me to an answer I gave to a similar question :
http://www.mapleprimes.com/questions/128128-Is-There-A-Way-To-Have-Filled-Or-Coloured

Using that idea we can do:

restart;
p:=proc(T) local p1,p2,A;
  p1:=plots:-polarplot([2-2*sin(theta),theta, theta = 0..T],thickness=3);
  A:=plottools:-getdata(p1)[-1];
  p2:=plots:-polygonplot(< <0|0>, A>,color=blue,style=patchnogrid);
  plots:-display(p1,p2)
end proc:

plots:-animate(p,[T],T=0..2*Pi);

In 1998 I wrote a procedure called filledplot. It took its basic idea from a procedure betweenplot written by Robert Israel.
I have updated it and here it is:

filledplot2.mw

You could combine two animations:
with(plots):
p1:=animate(plot,[ [[0,0],[2-2*sin(theta),theta]],coords=polar,axiscoordinates=polar ],theta=0.. 2*Pi,frames=250,scaling=constrained,trace=249):
#p2:=animatecurve([2-2*sin(theta),theta,theta=0.. 2*Pi],coords=polar,axiscoordinates=polar,frames=250, numpoints=300, filled=true,scaling=constrained):
p2:=animatecurve([2-2*sin(theta),theta,theta=0.. 2*Pi],coords=polar,axiscoordinates=polar,frames=250, numpoints=300,scaling=constrained):
display(p1,p2);