This works:

a:=t^2;

f:=unapply(a, t);

                                     a := t^2

                                  f := t -> t^2

Example:

f(3);

                 9

Definition of the function and its plotting for  alpha=2:

f := (x, y, alpha)->piecewise((-(1/2)*Pi <= x and x < 0 or 0 < x and x <= (1/2)*Pi) and -abs(x) <= y and y <= abs(x), (x^2+y^2)^alpha*arcsin(y/x),  x = 0 and y = 0, 0,  undefined);

plot3d(f(x, y, 2), x = -3 .. 3, y = -3 .. 3, numpoints = 25000, style = surface, axes = normal, orientation = [30, 70]);

Changing  the number  alpha  you will receive a variety of functions and their graphs.

restart;

with(plots):

B := 1.5:

u2 := r*cos(th): v2 := r*sin(th):

w1 := u1+I*v1: w2 := u2+I*v2:

z1 := evalc(w1^2): z2 := evalc(w2^2):

x1 := evalc(Re(z1)): x2 := evalc(Re(z2)):

y1 := evalc(Im(z1)): y2 := evalc(Im(z2)):

f1 := proc (theta) options operator, arrow; plot3d([-x1, -y1, v1], u1 = -6 .. 1, v1 = -B .. B, grid = [50, 50], orientation = [theta, 80], color = u1) end proc:

f2 := proc (theta) options operator, arrow; plot3d([-x2, -y2, v2], r = 1 .. 1, th = -Pi .. Pi, grid = [50, 50], orientation = [theta, 80], color = black, thickness=2) end proc:

display(seq(display(f1(10*k), f2(10*k)), k = -17 .. 18), insequence = true, axes = box, view = [-1 .. 1, -1 .. 1, -B .. B]);

eq:=diff(x(t),t,t) - 3*diff(x(t),t) + 2*x(t) = 2*t -3;

inc:=x(0)=2, D(x)(0)=4;

dsolve({eq, inc});

In radians Maple calculates immediately:

cos(Pi/6),  tan(Pi/4);

                                        1/2*3^(1/2),  1

If the angles are given in degrees, it is necessary to convert to radians:

cos(convert(30*degrees, radians)), tan(convert(45*degrees, radians));

                                        1/2*3^(1/2),  1

A, B, C := <x , 2 , -1>, <0, y , 1>, <3 , 4 , -2>:

AB, BC, CA := B-A, C-B, A-C:

d := u->LinearAlgebra[Norm](u,2,conjugate=false):

p := (u,v)->LinearAlgebra[DotProduct](u,v,conjugate=false):

Sol := [evalf(RealDomain[solve]({p(CA,AB)=0,d(CA)=d(AB)})),

evalf(RealDomain[solve]({p(AB,BC)=0,d(AB)=d(BC)})),

evalf(RealDomain[solve]({p(BC,CA)=0,d(BC)=d(CA)}))];

V:=map2(eval, [A,B,C], Sol):

T:=seq(plottools[polygon](map(convert,V[i],list)), i=1..4):

plots[display](T, axes=normal, scaling=constrained, orientation=[50,70]);

Edit: at first I did not notice that only needs to be AB = CA. Therefore, from the solutions should be left only the first two.

Example if  2 circles in the same plane:

eq1:=<0.8*cos(t),0.8*sin(t),0>:

eq2:=<2*cos(t),2*sin(t),0>:

C1:=plots[spacecurve](eq1, t=0..2*Pi, color=blue, thickness=3):

C2:=plots[spacecurve](eq2, t=0..2*Pi, color=blue, thickness=3):

S:=plot3d(s*eq1+(1-s)*eq2, t=0..2*Pi, s=0..1, style=surface, color=yellow):

plots[display](map(plottools[rotate],map(plottools[translate],[C1,C2,S], 3,4,1),Pi/3,0,Pi/6), axes=normal, orientation=[-15,60]);

 This method is suitable not only for the concentric circles and even if they lie in different planes:

eq1:=<0.8*cos(t)+0.5,0.8*sin(t),1>:

eq2:=<2*cos(t),2*sin(t),0>:

C1:=plots[spacecurve](eq1, t=0..2*Pi, color=blue, thickness=3):

C2:=plots[spacecurve](eq2, t=0..2*Pi, color=blue, thickness=3):

S:=plot3d(s*eq1+(1-s)*eq2, t=0..2*Pi, s=0..1, style=surface, color=yellow):

plots[display](map(plottools[rotate],map(plottools[translate],[C1,C2,S], 3,4,1),Pi/3,0,Pi/6), axes=normal, orientation=[120,45], scaling=constrained);

I think that to solve your problem it is convenient to use  Optimization  package:

restart;

de1 := diff(V(t), t) = V(t)-(1/3)*V(t)^3-W(t)+Ie:

de2 := diff(W(t), t) = 0.8e-1*(V(t)+.7-.8*W(t)):

sys:={de1, de2}:

sol:=dsolve({op(eval(sys,Ie=0.5)), V(0)=0, W(0)=1}, numeric, output=listprocedure);

Optimization[NLPSolve](rhs(sol[2]), 0..200, maximize, method=branchandbound );

Optimization[NLPSolve](rhs(sol[2]), 0..200);

beta := 1: lambda := 5: nu := 1: Delta := 1: 

sys := {C(0) = 0, In(0) = .2, diff(C(u), u) = beta*In(u)*s(u)-C(u)/nu, diff(In(u), u) = C(u)/nu-In(u)/lambda+Delta*(diff(diff(In(u), u), u)), diff(s(u), u) = -beta*In(u)*s(u), s(0) = .8, (D(In))(0) = .2}: 

dsys := dsolve(sys, numeric, [In(u), s(u), C(u)]): 

plots[odeplot](dsys, [[u,In(u)], [u,s(u)], [u,C(u)]], u = 0 .. 0.5, thickness=2, axes = boxed, legend=[In(u),s(u),C(u)]);

Because part of the plot is not included in the domain of the function  x, we extend the plot by  line segments using piecewise  command.

restart;

Eq1 := x = phi*(theta[m]-1/x)/theta[m]+(1-phi)*(1/2):

Sol := solve(Eq1, x):

x := unapply(Sol[1], phi):

p1 := plot(piecewise(`and`(theta[m] > .9, theta[m] < 1), 5*theta[m]-4.5, theta[m] > 1, x(0)), theta[m] = .9 .. 14, color = red, legend = "a=0 (extreme bimodal)"):

a := eval(x(.1), theta[m] = 1.55):

p2 := plot(piecewise(`and`(theta[m] > .9, theta[m] < 1.55), a*(theta[m]-.9)/(.65), theta[m] > 1.55, x(.1)), theta[m] = 1 .. 14, color = green, linestyle = 3, legend = "a=0.1"):

b := eval(x(.5), theta[m] = 3.56):

p3 := plot(piecewise(`and`(theta[m] > 3.55, theta[m] < 3.56), b*(theta[m]-3.55)/(0.1e-1), theta[m] > 3.56, x(.5)), theta[m] = 3.55 .. 14, color = blue, linestyle = 6, legend = "a=0.5"):

p4 := plot(piecewise(`and`(theta[m] > 3.9, theta[m] < 4), 5*theta[m]-20, theta[m] > 3.9, x(1)), theta[m] = 3.9 .. 14, color = black, linestyle = 2, legend = "a=1 (uniform)"):

P1 := plot([map(proc (t) options operator, arrow; [t, eval(x(0), theta[m] = t)] end proc, [1, 2, 8, 12])], style = point, color = red, symbol = cross, symbolsize = 17):

P2 := plot([map(proc (t) options operator, arrow; [t, eval(x(.1), theta[m] = t)] end proc, [1.55, 2.5, 3.5, 8, 12])], style = point, color = green, symbol = diagonalcross, symbolsize = 17):

P3 := plot([map(proc (t) options operator, arrow; [t, eval(x(.5), theta[m] = t)] end proc, [3.56, 5.5, 6, 8, 12])], style = point, color = blue, symbol = asterisk, symbolsize = 17):

P4 := plot([map(proc (t) options operator, arrow; [t, eval(x(1), theta[m] = t)] end proc, [4, 5, 8, 12])], style = point, color = black, symbol = box, symbolsize = 15):

plots[display]({P1, P2, P3, P4, p1, p2, p3, p4}, axes = boxed, view = [0 .. 14, 0 .. 1]);

p := [[565,575],[685,595],[700,580],[770,610]];

P:=plots[pointplot](p, color=red, thickness=2, connect=true):

T:=plots[textplot]([seq([op(p[i]),i],i=1..4)],align=right, font=[TIMES,ROMAN,18]):

plots[display](P,T);

Consider the function defined by the right-hand side of your equation and its typical plot (for b = 1)

f:=(z,b)->(arctan(b^2/(z*sqrt(b^2+b^2+z^2)))+b^2*z*1/(z^2+b^2)+1/(z^2+b^2))/sqrt(b^2+b^2+z^2))/(2*Pi):

plot(f(z,1), z=-3..3, thickness=2, discont);

The plot clearly shows the following properties of the function that are easy to prove rigorously:

1) The function is odd

2) The domain of the function is the all real numbers excepting  0

3) If  b<>0  then limit(f(z,b), z=0, right)=1/4  and  limit(f(z,b), z=infinity)=0

4) For  z>0  the function  f(z,b)>0  and it is strictly decreasing. This follows from the value of the derivative:

simplify(diff(f(z,b), z));

On the basis of the mentioned properties is easy to write a procedure that solves the equation for any a and b:

Sol:=proc(a,b)

if b=0 then if a=0 then return RealRange(-infinity,Open(0)), RealRange(Open(0),infinity) else return `No solutions` fi; fi;

if abs(a)>=1/4 or a=0 then return `No solutions` else

fsolve(f(z,b)=a, z) fi;

end proc:

Examples of use:

Sol(0.1, 2);

Sol(1, 2);

Sol(-0.1, 2);

Sol(0, 2);

Sol(1, 0);

Sol(0, 0);

For the beautiful plot use in addition to the option  discont=true  also  the option  symbol=solidcircle  and some other options:

plot([floor(x), seq([i, t, t = min(0, i) .. max(0, i)], i = -3 .. 3)], x = -3 .. 3.9, thickness = [3, 1$7], color = red, linestyle = [1, 2$7], discont = true, symbol = solidcircle, symbolsize = 15);

N := 10; h := 1/N;

for i from 0 to N do x[i] := i*h; p[i] := evalf(cos(x[i])) end do:

f1 := [seq([x[i], p[i]], i = 0 .. N)];

plot([cos(x), f1], x = 0 .. 1, y = 0 .. 1, style=[line,point], color=[red,blue],symbolsize=15);

Instead of  from i = 1 to N do   you should write  

for i from 1 to N do  

or  

for i to N do  

(initial value  i  equals 1 by default)