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)

For the rigorous proof of the absence of discontinuities in your example, you can use the identity

If the functions  f(x)  and  g(x)  are continuous,  this identity implies that  max(f(x), g(x))  is also continuous, because the function  abs  is continuous.

Obviously, the system  sys:= [eqn1-eqn2 = 0, eqn2-eqn3 = 0]  has the trivial solution  x=y=z=t  for any real  t . If I understand your goal, then you are trying to prove that this solution is the unique solution of the system. But this is not true, it is easy to see from the plot:

z := 1: 

sys:= [eqn1-eqn2 = 0, eqn2-eqn3 = 0]:

plots[implicitplot](sys, x = 0 .. 4, y = 0 .. 4, color = [red, blue], thickness = 2, gridrefine = 3);

My hypothesis  -  in typical case (excluding some individual values  t ) for each value of  z=t  the system  sys  has 4 solutions:[x=t, y=f(t), z=t],  [x=f(t), y=t, z=t],  [x=t, y=t, z=t], [x=g(t), y=g(t), z=g(t)]  where  f  and  g  are some functions.