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));

I know only the way by  textplot  command:

test := proc(x)

local A, B;

A:=plots[textplot]([0,0, `x is equal to:  `], color=green, font=[TIMES, ROMAN, 18], align=left);

B:=plots[textplot]([0,0, 2], color=black, font=[TIMES, ROMAN, 18], align=right);

plots[display](A, B, axes=none);

end proc:

Example:

test(2);

test := proc (x)

cat(`x is equal to:  `, x);

end proc:

Example:

test(2);

                x is equal to:  2

A:=proc(M)

Matrix(M, M, {(1,1)=1/2, (1,2)=1/2/sqrt(2), (2,1)=-1/8/sqrt(2), seq((i,1)=-1/2/sqrt(2)/i/(i-2), i=3..M), seq((i,i-1)=-1/4/(i-2), i=3..M), seq((i,i+1)=1/4/i, i=2..M-1)});

end proc:

Examples:

A(3), A(4), A(5);

If the parameter values are specified​​, the equation maybe easily solved numerically by finding all the roots. To simplify the typing I changed the notations of variables and parameters. At the output all are unchanged.

Example:

x:=alpha: t0:=mu0: t1:=mu1: s:=sigma:

eq:=60*t0*x^2*s^2-60*x*s^2*t0^2+60*t1*s^2*x^2-60*s^2*t1^2*x-15*t0*x^2-15*t1^2*x^4+30*t1^2*x^3-30*t1^3*x^2+15*t1^4*x-40*x^3*s^2+30*x^3*t0^2-30*t0^3*x^2+240*x*s^4+15*x*t0^4+6*x^5+20*s^2-60*t1*s^2+60*t1^2*s^2-120*s^4+15*t1-30*t1^2+30*t1^3-15*t1^4-3=0; 

fsolve(eval(eq, {mu0=1, mu1=2, sigma=3}), alpha, complex);