restart:  with(linalg):
S:=matrix(3,3,[s11,s12,s13,s21,s22,s23,s31,s32,s33]):
X:=matrix(3,1,[x1,x2,x3]):
A:=matrix(3,3, (i,j)->S[i,j]/X[j,1]);

Xm:=evalm(transpose(h)&*transpose(X));
S/~Xm;  # memberwise division: S / Xm

Xd:=diag(x1,x2,x3);  # make X the diagonal of a square matrix
A:=evalm(S&*inverse(Xd));
 

Example:

P_X:=Vector(10, symbol=v);

seq(P_X[i]<1, i=1..10);

Addition for beginners. The command  seq  is one of the most powerful and frequently used commands in Maple. It can often substitute for loops, both simple and nested (in the latter case, nested sequences are used). Here are 2 more examples of using this command:

seq(ithprime(n), n=1..50);   # First 50 prime numbers

n:=50:
seq(seq(seq(`if`(a^2+b^2=c^2, a^`2`+b^`2`=c^`2`, NULL), b=a+1..n-1), a=1..n-2), c=3..n);   # All Pythagorean triangles with hypotenuse not exceeding 50

Edit.

I do not think that this sum can be found in closed form for any  r . But for integer positive  r  the solution can be written as a procedure:

restart;
S:=r->simplify(sum((a+k*d)^r, k=1..n));

Examples of use:

S(2);
S(3);
S(6);

See  wiki  for details.

Example of calculation (norm in the space  C[-1,1]):

restart;
f:=x->x^3-x:
maximize(abs(f(x)), x=-1..1);
                                                        2*sqrt(3)*(1/9)

NumberTheory:-Totient(n)  returns the number of positive integers coprime to  n  and not greater than  n . See help for details.
 

Addition. Another way:

restart;
with(Physics):
phi:=n->add(KroneckerDelta[igcd(n,k),1], k=1..abs(n)-1);
NumberTheory:-Totient(14);
phi(14);
 

Edit.
 

See  here  about this.

Of cause you can just write  3+1/7 .

Here is an example of the 2 ways to calculate the value of the derivative of a function at a point. Your function can be specified as a procedure or as some expression. In the first case, it is convenient to use the differential operator  D , in the second case  eval  command:

restart;
f:=(x,y)->x^2+x*y^3:  # The function is given as a procedure
g:=x^2+x*y^3:   # The function is given as an expression
D[1](f)(0.5, 1);   # The value of the partial derivative with respect to x at the point
eval(diff(g, x), [x=0.5, y=1]);  # The same
                                                     2.0
                                                     2.0

Below the example of the same error that you have. You can not differentiate by number:

x:=0.5:  y:=1:
diff(g, x);
   Error, invalid input: diff received .5, which is not valid for its 2nd argument
 

I just reduced the domain of the function:

with(plots):

animate(plot3d,[sin(Pi*(t-x))*sin(Pi*y)*Heaviside(t-x)-sin(Pi*(t+x))*sin(Pi*y)*Heaviside(t+x),x=0..4,y=piecewise(x>=0 and x<2,1/2-x/4,0)..piecewise(x>=0 and x<2,1/2+x/4,1)],t=-2*Pi..2*Pi, frames=90);
 

You can do something similar to this simple example:

restart;
for y from 1 to -1 by -0.1 do
sol:=fsolve(exp(x)=y);
if type(sol,numeric)=false then break fi;
od;
y;

As soon as y = 0 the calculation is interrupted.

restart;
f := (x,y)->sin(x)*cos(y):
plot3d(f, 0..2*Pi, 0..2*Pi);
plots:-contourplot(f, 0..2*Pi, 0..2*Pi);

PS. If you want contours with numeric labels, then see this post

restart;
sol := dsolve({(diff(r*(diff(u(r), r)), r))/r = 0, u(R0) = U2, u(h) = U1}, u(r));
plot(eval(rhs(sol), [h = 1, R0 = 5, U1 = 2, U2 = 3]), r = 1 .. 5);

It is much easier to calculate the length of this arc of the ellipse (half of the total ellipse) using its parametric equations:

eqn := (1/32)*(x-16)^2+(1/2025)*(y+0)^2 = 1;
a:=sqrt(32):  b:=sqrt(2025):
int(sqrt(diff(16+a*cos(t),t)^2+diff(b*sin(t),t)^2), t=0..Pi);
evalf(%);
plots[implicitplot]([(1/32)*(x-16)^2+(1/2025)*(y+0)^2 = 1], x = 0 .. 32, y = 0 .. 60, scaling = constrained);

To plot arrows, use  plots:-arrow  command, and for text use  plots:-textplot command. See help on these commans for details.

Example of the plotting:

with(plots):
P:=plot([seq(5-a/20*(x+5)^2, a=[1,2,3])], x=-5..3, -2..7, color=[red,blue,green], size=[700,300]):
A1:=arrow([-2,6],[-0.5,-1.85], width=[0.01,relative=false], head_width=[0.07,relative=false], head_length=[0.2,relative=false]):
A2:=arrow([-1.3,6],[-0.5,-1.9], width=[0.01,relative=false], head_width=[0.07,relative=false], head_length=[0.2,relative=false]):
A3:=arrow([-0.6,6],[-0.5,-1.65], width=[0.01,relative=false], head_width=[0.07,relative=false], head_length=[0.2,relative=false]):
T:=textplot([[-0.6,6.4,a=1],[-1.3,6.4,a=2],[-2,6.4,a=3]], font=[times,14]):
display(P,A1,A2,A3,T);

Your integral contains a large number of parameters. Integral is easy to calculate if you specify the values of these parameters. Also note that you should not use a variable and an indexed variable  with the same name in the same worksheet (you have d and d[i], x and x[i]).

Example:

restart;
eval(sum(d[i]*ln(x[i]/theta), i = 1 .. 3)/(d1*(sum(ln(x[i]/theta), i = 1 .. 3))^v),[v=2,d1=6,d[1]=1,d[2]=2,d[3]=3,x[1]=1.5,x[2]=2.3,x[3]=3.7]);
int(%, theta=0..2);
                          

Try this version. In classic GUI of Maple 2017.3 it works:

restart;
with(plots):
with(plottools):
w:=1.5:
lup := line([w,0.9], [w,5.1], color=red, thickness=4):
plt:=plots[implicitplot]((x-1)^4/24+(y-3)^4/10=2, x=w..7, y=0..12,color=red, thickness=4):
op([1,1],plt);
convert(%,listlist):
P:=polygon(%, color=gold):
plots[display]([lup,plt],scaling=constrained);
plots[display]([lup,P,plt],color=gold,scaling=constrained);