evalf(Int(sqrt(x^3+1), x=1..2));

2.129861293

I think I understand what you need. Use a double loop as follows:

L:=[0,1,2]:

M:=[]:

for i to nops(L) do

for j to nops(L) do

M:=[op(M), [L[i], L[j]]]:

od: od:

M;

If your polynomial is a quadratic form, then use the Sylvester criterion. For an arbitrary polynomial  the appropriate command in Maple is unknown to me.

Repetitive elements can easily remove from the list, making the set from the list.

L:=[0, 0, 1, 1, 2, 2]:

M:={op(L)}:

combinat[choose](M,2);

{{0,1},{0,2},{1,2}}

a:=Matrix([[4,3,4],[7,7,3],[7,3,1]]):

b:=Matrix([[6,3,5],[9,8,0],[8,8,1]]):

L:=[a,b,b,a]:

add(LinearAlgebra[DiagonalMatrix]([0$(2*i-2),L[i],0$(8-2*i)]),i=1..nops(L));

r:=25:

Sol:=dsolve({diff(x(t),t)=-10*x(t)+10*y(t), diff(y(t),t)=r*x(t)-y(t)-x(t)*z(t), diff(z(t),t)=(-8/3)*z(t)+x(t)*y(t), x(0)=-1, y(0)=-1, z(0)=1}, numeric):

plots[odeplot](Sol,[x(t), y(t), z(t)], 0..30, axes=normal, numpoints=10000);

It turns out a very interesting picture:

See  plots[matrixplot]  command.

f := (x, y) -> ((2*x^2+y^2*(x^2*(-2+sqrt(1/x^4+4/(x^2*y^2)+4/y^4+4))-1))/(4*x^2-2*y^2+4))^(1/2) ;

g := unapply(diff(f(x, y), x), x,y);

plot([g(.1, y), g(1, y), g(3, y)], y = 0 .. 10);

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

interface(rtablesize=100):

B:=LinearAlgebra[DiagonalMatrix]([A, -A, A]);

Instead of L in fsolve command after comma write L=0..infinity

You have a rather complicated transcendental equation. Solve command does not solve such equations. Replace solve command by fsolve command.

A:=plot(x, x=0..1):

B:=plot(x^2, x=0..1):

C:=plot(x^5, x=0..1):

plots[display](array([A, B, C]), scaling=constrained, thickness=2);

The function

x->2*a/Pi*arcsin(sin(2*Pi*x/b))

gives a triangular signal between -a and a with the period b .

An example:

a:=1: b:=4:

plot(2*a/Pi*arcsin(sin(2*Pi*x/b)), x=0..8, thickness=2, scaling=constrained);

Dr:=[[1,5],[2,6],[3,4],[4,5],[5,7],[6,8],[7,4],[8,6],[9,5],[10,7]]:

Y:=proc(i) Dr[i,2] end; ty:=proc(i) Dr[i,1] end;

P[X](omega)=(1/2)*{(Sum('Y'(j)*cos(omega*('ty'(j)-1/(2*omega*tan((Sum(sin(2*omega*'ty'(j)),j=1..N[0]))/(Sum(cos(2*omega*'ty'(j)),j=1..N[0]))))))^2,j=1..N[0]))/(Sum(cos(omega*('ty'(j)-1/(2*omega*tan((Sum(sin(2*omega*'ty'(j)),j = 1..N[0]))/(Sum(cos(2*omega*'ty'(j)),j=1..N[0]))))))^2,j= 1..N[0]))+(Sum('Y'(j)*sin(omega*('ty'(j)-1/(2*omega*tan((Sum(sin(2*omega*'ty'(j)),j=1..N[0]))/(Sum(cos(2*omega*'ty'(j)),j=1..N[0]))))))^2,j=1..N[0]))/(Sum(sin(omega*('ty'(j)-1/(2*omega*tan((Sum(sin(2*omega*'ty'(j)),j=1..N[0]))/(Sum(cos(2*omega*'ty'(j)),j =1 .. N[0]))))))^2,j=1..N[0]))};

Walk1 := proc(n)

local pick, i,j,edge,step,L,a, A, B, C, E;   

pick := rand(0..1);    (i,j) := (0,0);  L:=[[0,0]];  edge := n/2;   

for step do a:=pick();       

if a=0 then i := i + 2*pick()-1; L:=[op(L), [i,j]]; fi;       

if a=1 then j := j + 2*pick()-1; L:=[op(L),[i,j]]; fi;       

if edge < abs(i) or edge < abs(j) then  print(step); break;       

end if; end do;

print(L);

A:=plot([seq([t, i, t=-n/2..n/2],i=-n/2..n/2), seq([i, t, t=-n/2..n/2],i=-n/2..n/2)], color=black);

B:=seq(plottools[disk](L[i],0.1, color=red), i=1..nops(L));

E:=seq(plottools[curve]([seq(L[k],k=1..i)], color=blue, thickness=5), i=1..nops(L));

C:=seq(plots[display](op([A, B[i], E[i]]))$10, i=1..nops(L));

print(plots[display](seq(C[i],i=1..10*nops(L)),insequence=true,scaling=constrained,view=[-n/2-1..n/2+1, -n/2-1..n/2+1]));

end proc:

The procedure returns the number of steps, a list of nodes traversed, and the animation of the random walk.

An example:

Walk1(4);