LengthOfIntersection procedure finds in the form of an explicit piecewise function the length intersection of 2 segments (real closed intervals) (see the original question). Now we can not only find this length for a given x, but also differentiate, integrate, etc.

LengthOfIntersection := proc (x::name, A::list(realcons), B::list(realcons))

local h, B1, d;

h := min(A[2]-A[1], B[2]-B[1]);

d := abs(B[2]-B[1]-A[2]+A[1]);

B1 := B+~x;

piecewise(convert(solve(B1[2] <= A[1], {x}), `and`), 0,

convert(solve({B1[2] < A[1]+h, A[1] < B1[2]}, {x}), `and`), B1[2]-A[1],

convert(solve({B1[2] <= A[1]+h+d, A[1]+h <= B1[2]}, {x}), `and`), h,

convert(solve({B1[1] < A[2], A[1]+h+d < B1[2]}, {x}), `and`), A[2]-B1[1],

convert(solve(A[2] <= B1[1], {x}), `and`), 0);

end proc:

Examples of use:

f := unapply(LengthOfIntersection(x, [1, 3], [0, 5]), x):

f(x);

f(2);

diff(f(x), x);

int(f(x), x = -4 .. 3);

plot(f, -5..5, color=red, thickness=3, scaling=constrained);

You can do 

for i to 4 do i^2 end do;

i;

i := 'i'; 

i;

Now  i  is symbol

whattype(i);

                                  symbol

Alternative option:

V:=[a5, a6, a7]:

map(t->[t,0], V):

seq(seq(seq([i,j,k], k=%[3]), j=%[2]), i=%[1]);

m:= [01100101, 01101100, 01100111, 01100001];

map(t->[seq(parse(convert(t,string)[i]), i=1..length(convert(t,string)))], m);

 or

map(t->parse~([StringTools[LengthSplit](convert(t, string),1)]), m);

Addition: If you want to preserve leading zeros, then the original numbers should be directly defined as  strings or symbols, for example:

m:= convert~([`01100101`, `01101100`, `01100111`, `01100001`], string);

map(t->[seq(parse(t[i]), i=1..length(t))], m);

Otherwise Maple automatically discards these zeros.

Example:

ListTools[Reverse]([1, 1, 0, 0, 1, 1, 1, 0]);

                                                 [0, 1, 1, 1, 0, 0, 1, 1]

Simple procedure  f  solves your problem:

f:=proc(x, A, B)

local B1;

B1:=[B[1]+x, B[2]+x];

`if`(A[2]<=B1[1] or B1[2]<=A[1], 0, min(A[2],B1[2]) - max(A[1],B1[1]));

end proc:

Examples of use:

f(0,[2,4],[0,1]);

f(3/2,[2,4],[0,1]);

f(3,[2,4],[0,1]);

f(4,[2,4],[0,1]);

f(0,[1,3],[0,5]);

plot('f'(x,[1,3],[0,5]), x=-5..5, color=red, thickness=3, scaling=constrained);

I don't know any Maple procedures which solve the similar problems.

1) LinearAlgebra[LinearSolve]  command finds the infinite number of solutions that depend on the three parameters:

restart;

with(LinearAlgebra):

 l1 := [1, 1, 1, 0, 0, 0, 0, 0, 0]:

l2 := [0, 0, 0, 1, 1, 1, 0, 0, 0]:

l3 := [0, 0, 0, 0, 0, 0, 1, 1, 1]:

l4 := [1, 0, 0, 1, 0, 0, 1, 0, 0]:

l5 := [0, 1, 0, 0, 1, 0, 0, 1, 0]:

l6 := [0, 0, 1, 0, 0, 1, 0, 0, 1]:

l7 := [0, 0, 1, 0, 1, 0, 1, 0, 0]:

A := Matrix([l1, l2, l3, l4, l5, l6, l7]):                                   

b := <15, 15, 15, 15, 15, 15, 15>:  # Vector specified by right hand side of the matrix equation  A.m=b

m:=LinearSolve(A, b);

2) The second problem has 24 solutions:

S:={$ 2..8}:  n:=0:  T:=table():

for i in S do

for j in S do

for k in S do

m0:=eval(m,{_t[6]=i, _t[8]=j , _t[9]=k});

if convert(m0, set)=S then n:=n+1; T[n]:=m0  fi;

od: od: od:

T:=convert(T, set);

nops(T);

 3) The third problem has no solutions:

S1:={$ 0..8}:

P:=combinat[permute](S1, 3): n:=0: T:=table():

for p in P do

m0:=eval(m,{_t[6]=p[1],_t[8]=p[2],_t[9]=p[3]});

if convert(m0,set)=S1 then n:=n+1; T[n]:=m0 fi;

od:

T:=convert(T,set);

                                                       T := {}

But if we replace  {$ 0..8}  by  {$ 1..9)  then again get  24  solutions:

S2:={$ 1..9}:

P:=combinat[permute](S2, 3):  n:=0:  T:=table():

for p in P do

m0:=eval(m, {_t[6]=p[1],_t[8]=p[2],_t[9]=p[3]});

if convert(m0, set)=S2  then n:=n+1; T[n]:=m0  fi;

od:

T:=convert(T, set):

nops(T);

                                                         24

Matrix_Equation.mws

In fact, you just need to plot a function by some points. Here is an example  of such a plotting. This is not a procedure, but the example of a general nature and can easily be re-written as a procedure:

restart;

f := unapply(x^2, x):  a := 1:  b := 3:  n := 10:

x := a+k*(b-a)/n:

L := [seq([x, f(x)], k = 0 .. n)]:

plot(L, color = red, thickness = 2);

If you type  I__K  you get just indexed symbol, to which you can assign any object, for example the integer  5 . If now you type  I[K]  this is not equal   I__K .  You get some strange nature of the object, you can not assign anything to it:

The reason is that in Maple  I  is not a symbol but the complex constant, namely, the imaginary unit. See

whattype(I);

I^2;

So, I[K]  never should be typed , because it does not make sense.

Usually a system of equation is specified as a set, because the order of equations is irrelevant.

L := [a, c, b, d, f, e]:   S := [b, c, f, e]:
`union`(op(map(convert, [L, S], set)));

                               {a, b, c, d, e, f}

Example:

interface(rtablesize=infinity):  # It lets to display matrices of arbitrary sizes
LinearAlgebra[RandomMatrix](15);

st := time():

for i to 1000 do for j to 1000 do

a[i, j] := abs(i-j+1)^0.3-abs(i-j)^0.3;

end do end do:

time()-st;

                                29.812

For example:

abs(5)^0.3-abs(4)^0.3;

                            0.104940030

I do not think that 29 seconds is too much time for this task.

My comp: 1.7GHz, 4.00Gb Ram

Procedure  DerPlot  solves all the problems. I added the plot of the derivative  f'(x) (yellow color), because the point  [a,f'(a)] (green color)  lies on the plot of the derivative  rather than on the plot  f(x) .

Formal parameters:  Expr  is an expression for  f(x) ,  a  is the point of tangency,  h (optional) is a numeric that specifies a range for  x  (x=a-h .. a+h) (by default  x=4) .

DerPlot := proc (Expr, a, h := 4)

local f, T;

uses plots;

f := unapply(Expr, x);

T := f(a)+(eval(diff(f(x), x), x = a))*(x-a);

plot([f(x), T, diff(f(x), x)], x = a-h .. a+h, color = [blue, red, yellow], thickness = 2, caption = typeset("Curve f(x) :  ", y = f(x), "\nTangent line:  ", y = T, "\n Red point :  ", [parse("a"), f(parse("a"))] = [a, f(a)], "\n Green point :  ", [parse("a"), `f'`(parse("a"))] = [a, eval(diff(f(x), x), x = a)]), captionfont = [times, roman, 16]);

pointplot([[a, f(a)], [a, eval(diff(f(x), x), x = a)]], color = [red, green], symbol = solidcircle, symbolsize = 16);

display(%, %%);

end proc:

Example of use:

DerPlot(sin(x), 2*Pi*(1/3));

Procedure  PartialSum calculates partial sums of the series   a1+a2+...+an+...

Formal parameters:  N  is the number of  terms, an is an expression for nth term of the series.

PartialSum:=proc(N, an)
local psum, f, k;

psum:=0;
f:=unapply(an, n);
for k from 1 to N do
psum:=psum+f(k);
od;
psum;

end proc:

Example of use:

PartialSum(10, 1/n);

                               7381/2520

Addition:  Of course,  PartialSum  procedure is simply a correction of  OP's attempt. In actual practice, you can simply write 

 add(an, n=1..N);  # for numeric  N

sum(an, n=1..N);  # for symbolic  N

Examples:

add(1/2^n, n=1..10);

sum(1/2^n, n=1..N);

eval(%, N=10);

Markiyan Hirnyk  In fact, you've got the original parameterization  x=9z^2+3z+367, y=6z^2 +z+244  and nothing new.

If  x, y  are integers, we will prove that there must be  z  is an integer (this is obvious)  or  z  is  an integer+1/3

Let   x, y  are integers. It is easy to get  3*y-x=9*z^2+365 . It follows that 3*z must be an integer. Let  3*z=n ,  z=n/3 , n is integer.  It is clear that  x  be  an integer for all  z=n/3 , n is integer. It remains to find conditions on n, y to be an integer.  To do this, we have to solve in integers  the equation  6*(n/3)^2+n/3=m :

isolve(6*(n/3)^2+n/3=m);

             {m = 6*_Z1^2-_Z1, n = -3*_Z1}, {m = 6*_Z1^2-5*_Z1+1, n = 1-3*_Z1}

From the second solution we have  z=n/3=-_Z1+1/3