1) Your code is very inefficient. To specify the sequence of the objects is better to use the  seq  command. Also I do not understand the meaning of the options  coords=polar, axiscoordinates=polar, as your circles given in Cartesian coordinates.

with(plots): with(plottools):

c := seq(circle([i/(1+i), 0], 1/(1+i), color = green), i=1..20):

d := seq(circle([1, 1/i], 1/i, color = blue), i=1..20);

display(c, d, view = [0 .. 2, -0.5 .. 2], scaling=constrained);

2) intersect command works only with finite sets. In order to find the point of intersection of the two lines, you can use solve  command.

Sol := solve({(x-1)^2+(y-1)^2 = 1, (x-2/3)^2+y^2 = 1/9}, explicit);

                        {x = 1, y = 0}, {x = 2/5, y = 1/5}

Point1 := disk([rhs(Sol[1, 1]), rhs(Sol[1, 2])], 0.01, color = red):

Point2 := disk([rhs(Sol[2, 1]), rhs(Sol[2, 2])], 0.01, color = red):

display(L, LL, Point1, Point2, view = [0 .. 1.5, -0.5 .. 0.5], scaling = constrained);

Edited. Fixed equation of  c .

Your equation has infinitely many solutions for any real  n>0:

RealDomain[solve](3^x/2^y = n, {x, y});

              {x = (ln(n)+ln(2)*y)/ln(3), y = y}

As  y  you can take any number, ie the set of solutions depends on a single parameter.

You can easily place this idea as a procedure. Procedure  Round  rounds each number float type in expression  expr   to  n  digits after decimal point.

Round:=proc(expr, n)

local F;

uses StringTools;

F:=x->parse(sprintf(Substitute("%.nf", "n", convert(n, string)), x));

subsindets(expr, float, F);

end proc;

Examples:

Round(Vector([1.235, 0.1237]), 2);

Round(1.156*x^2-12*x+0.533, 1);

See acer's answer here

a:=0.0000000000000000000000213123123:

parse(sprintf("%.3f", a));

                                        0.000

Plots := proc (L, C)  # L - the list of exponents, C - the list of colors

local k, f, M;

f := (x, n)->1+2*(1+x)^n;

M := mul(f(5, k), k = L)^(1/nops(L));

print(plots:-display(Array([seq(plot(f(x, L[k]), x = 0 .. 5, 0 .. M, color = C[k], caption = typeset("The exponent   ", n = L[k])), k = 1 .. nops(L))])));

plot([seq(f(x, k), k = L)], x = 0 .. 5, 0 .. M, color = C, legend = [seq(n = k, k = L)]);

end proc;

Example:

Plots([0, 1, 2, 3], [yellow, green, blue, red]);

1) To emphasize the symmetries of the graphs, the center of the composition take the point (-1,1).

2) For greater clarity, the range along y-axis is reduced.

3) In order to estimate "by eye" the slopes of the graphs, selected the same scales on axes.

plot([seq(1+2*(1+x)^n, n = 1 .. 4)], x = -3 .. 1, -2 .. 4, color = [red, blue, yellow, green], thickness = 2, legend = [n = 1, n = 2, n = 3, n = 4], scaling = constrained);

restart;

a:=Vector([2,3,4,5]): l:=[]:

for k to ArrayNumElems(a) do

if a[k]>=3 then l:=[op(l), k] fi:

od:

l;

                        [2, 3, 4]

Try first all equations lead to a uniform form.

Example:

A := y - x = 0:

k := lcoeff(lhs(A), [x, y]);

A/k;

                k := -1                

               x - y = 0

Carl, you are right. I was inattentive. New procedure   SpiralMatrix1  solves the original problem:

restart;

SpiralMatrix1:=proc(A)

local n, A1, SpiralMatrix, a;

uses LinearAlgebra;

n:=RowDimension(A);

SpiralMatrix:=proc(A::Matrix)

local n, Turn, L, C, k;

uses LinearAlgebra;

n:=RowDimension(A);

Turn:=proc(B::Matrix)

local n;

uses LinearAlgebra;

n:=RowDimension(B);

[[seq(B[1,j],j=1..n),seq(B[i,n],i=2..n),seq(B[n,j],j=n-1..1,-1),seq(B[i,1],i=n-1..2,-1)], Matrix([seq([seq(B[i,j],j=2..n-1)],i=2..n-1)])];

end proc;

L:=[]; C:=A;

for k to floor(n/2) do

L:=[op(L),op(Turn(C)[1])]; C:=Turn(C)[2];

od;

if type(n,even) then Matrix(n,n,L) else Matrix(n,n,[op(L),C[1,1]]) fi;

end proc; 

A1:=SpiralMatrix(Matrix(n,n,symbol=a));

assign(seq(seq(A1[i,j]=A[i,j],j=1..n),i=1..n));

Matrix([seq([seq(a[i,j],j=1..n)],i=1..n)]);

end proc;

Example:

A:=Matrix([seq([i $ 5], i=1..5)]);

SpiralMatrix1(A);

SpiralMatrix:=proc(A::Matrix)

local n, Turn, L, C, k;

uses LinearAlgebra;

n:=RowDimension(A);

Turn:=proc(B::Matrix)

local n;

uses LinearAlgebra;

n:=RowDimension(B);

[[seq(B[1,j],j=1..n),seq(B[i,n],i=2..n),seq(B[n,j],j=n-1..1,-1),seq(B[i,1],i=n-1..2,-1)], Matrix([seq([seq(B[i,j],j=2..n-1)],i=2..n-1)])];

end proc;

L:=[]; C:=A;

for k to floor(n/2) do

L:=[op(L),op(Turn(C)[1])]; C:=Turn(C)[2];

od;

if type(n,even) then Matrix(n,n,L) else Matrix(n,n,[op(L),C[1,1]]) fi;

end proc;

Example:

A:=Matrix([seq([i $ 5], i=1..5)]);

SpiralMatrix(A);

restart;

f(x) = 2.25*f(x-1)-0.5*f(x-2);

convert(%, rational);

f := unapply(rsolve({%, f(1) = 1/3, f(2) = 1/12}, f(x)), x);

plots[display](<plot(f, 1 .. 40, color = red) | plots[logplot](f, 1 .. 40, color = red)>);

Left - the usual plot, right - logplot.

Here is the simplest code, generating the first  N  rows of Pascal's Triangle:

PascalTriangle:=proc(N)

local n;

for n from 0 to N-1 do print(seq(binomial(n,k), k=0..n)) od;

end: 

Example:

PascalTriangle(10);

Line:=plot(x^2-2*x+2, x=-0.5..2.5, -0.5..3.5, thickness=2):

Point:=plottools[disk]([1,1],0.04, color=blue):

plots[display](Line, Point, scaling=constrained);

For such a simple task we don't need in Maple (and it will not help if, for example, instead of 1000 we take 10 ^ 10). Instead the range 1 .. 1000 consider the range 0 .. 999, adding zeros in front if the number of digits is less than 3:  000, 001, 002, etc. Total used 3 * 1000 = 3000 digits and because of the equality of all digits  the answer will be 3000/10 = 300

The procedure  MultiAdd  solves your problem:

MultiAdd:=proc(f::procedure, L::list(range))

local M, T, m, S;

uses combinat;

M:=map(i->[$ i], L);

T:=cartprod(M);  S:=0;

while not T[finished] do

m:=T[nextvalue]();  S:=S+f(op(m))

end do;

S;

end proc:

Example:

MultiAdd((i,j,k)->1/(i*j*k), [1..5, 2..8, 3..7]);

                              3360747/784000