You can use  RealRange command . For example:

is(2 in RealRange(0, Open(2)));

is(1 in RealRange(0, Open(2)));

is(0 in RealRange(0, Open(2)));

                              false

                              true

                              true

When solving equations, you can just use the corresponding inequalities. Because in your equation 3 unknowns, two of them must be specified. For example:

restart;

y:=2:  z:=0.3:

solve({y=1.048 + 1.02*x + 6.118*(z-4.041*x^2) + 16.22*(z^2) +6.241* (x*z),

x>=0.001, x<=0.7}, x);

                                                         {x =0 .3718638796}

a := x->piecewise(x^2 <= 3, x^2, x^1.5):

a(1);  a(2);

                  1

        2.828427125

Addition: one more way

a := x->`if`(x^2 <= 3, x^2, x^1.5):

convert(exp(x), Sum, dummy=n);

This is quite easy. For example:

V:=[V[k] $ k=1..n]: N:=[$ 1..n]:

plot(N, V, thickness=2);

Even for a specified  b the  solution is expressed only in terms of integrals:

for b from -5 to 5 do

P[b](t):=t*ln(t)^b:

dsolve(diff(Q(t), t)^2 = diff(P[b](t), t), Q(t)):

od;

1st column of  matrix  A  by  3rd column of matrix  B :

A := Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]]);
B := Matrix([[9, 8, 7], [6, 5, 4], [3, 2, 1]]);

A[.., 1] .~ B[.., 3];  whattype(%);

Edited: the extra bracket was deleted. 

restart;

with(LinearAlgebra):

A := `<|>`(`<,>`(1, 2,3), `<,>`(2, 3, 0), `<,>`(2, 0, 0));

v, EigenVector1:= Eigenvectors(A):

NewMatrix3 := Matrix([[x1,x2,x3], [x4,x5,x6], [x7,x8,x9]]):

simplify(solve({seq(seq((NewMatrix3 . EigenVector1(1..-1,i))[j]=(v[i]*EigenVector1(1..-1,i))[j], j=1..3), i=1..3)}, {seq(x||i, i=1..9)}));

assign(%);

NewMatrix3;

 Edited:  restart  command was added.

L := [3, -2, 5, 1, 10, 7]:

[seq(add(L[i], i = 1 .. k), k = 1 .. nops(L))];

                   3, 1, 6, 7, 17, 24

Addition:  if  L  is a long list, then more effective variant is

L := [3, -2, 5, 1, 10, 7]:  S[1] := L[1]:  n := nops(L):

for i from 2 to n do  S[i] := S[i-1]+L[i]  end do:

S := [seq(S[i], i = 1 .. n)];

I just built your plot by points (a total of 41 points were used):

f:=phi->Int(309.9615685*exp(3.936759889*phi-3.936759889*`&varphi;`)*Heaviside(-phi+`&varphi;`)*exp(-1549.807842*`&varphi;`^2)/(Pi*sqrt(-15.74703956*phi+15.74703956*`&varphi;`)), `&varphi;` = -infinity .. infinity);

X:=[seq(x, x=-1..1, 0.05)]:

Y:=[seq(f(X[i]), i=1..nops(X))]:

plot(X, Y, thickness=2);

assume(L > 0);
assume(K > 0);
assume(mu > 0);
assume(a<0);
assume(t > 0);
convert(cosh((1/2)*t*sqrt(a)/L), sincos);
subs(a=mu^2-4*L*K, %);

Code is edited.

You can manually make these axes using  plottools  package like this  http://www.mapleprimes.com/questions/200601-Labels-Above-Axis#answer202152

The first integral is calculated incorrectly (a bug). If test1 to split into two integrals, then all is well:

evalf(Int(cot(2*t-1-I), t=0..0.5))+evalf(Int(cot(2*t-1-I), t=0.5..1));

                                         0.+1.115976093*I

You can generate a random Hermitian matrix as follows:

LinearAlgebra[RandomMatrix](4, generator=rand(-9..9)+I*rand(-9..9), outputoptions=[shape=hermitian]);

A:=x/2+y/2:

subs(a=``, normal(A/a));