There is no triangle with sides 136, 170, 174, in which the coordinates of all vertices are integers or rational numbers. It is well known that the area of any such triangle is rational. But the area of your triangle is irrational:

a := 136: b := 170: c := 174:

p:=(a+b+c)/2:

sqrt(p*(p-a)*(p-b)*(p-c));

                                        240*2002^(1/2)

Tools -> Options -> Precision -> Round calculation ... -> Apply to Session or Apply Globally

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