The method you are referring to is based on the usual shift and has a significant drawback - the thickness of such a surface in different places will be different. This is clearly seen in 2D-examples below. In the first plot, the second curve is obtained by shifting the parabola  y=x^2  by 0.1 up. In the second graph, the second curve passes through the ends of the all normals with the length 0.1 at the points of the parabola  y=x^2 :

restart;
plot([x^2, x^2+0.1], x=-1.5..1.5, color=black, scaling=constrained);
A:=plot(x^2, x=-1.5..1.5, color=black, scaling=constrained):
B:=plot([x-2*x/sqrt(1+4*x^2)*0.1, x^2+1/sqrt(1+4*x^2)*0.1, x=-1.5..1.5], color=black, scaling=constrained):
plots:-display(A, B);

Below is the same method with normals applied to the surface defined by the explicit equation  z=x^2-y^2. For clarity, both sides of the surface and its end face are made in different colors:

f:=(x,y)->x^2-y^2:
A:=plot3d(f(x,y), x=-1..1, y=-1..1, color=pink, style=surface):
Dx:=diff(f(x,y),x): Dy:=diff(f(x,y),y):
d:=sqrt(Dx^2+Dy^2+1): h:=0.1:
B:=plot3d([x-Dx/d*h,y-Dy/d*h,f(x,y)+h/d], x=-1..1, y=-1..1, color=blue, style=surface):
C1:=plot3d([1-2/sqrt(2^2+Dy^2+1)*h*s,y-Dy/sqrt(2^2+Dy^2+1)*h*s,f(1,y)+h*s/sqrt(2^2+Dy^2+1)], s=0..1, y=-1..1, color=cyan, style=surface):
C2:=plot3d([-1+2/sqrt(2^2+Dy^2+1)*h*s,y-Dy/sqrt(2^2+Dy^2+1)*h*s,f(-1,y)+h*s/sqrt(2^2+Dy^2+1)], s=0..1, y=-1..1, color=cyan, style=surface):
C3:=plot3d([x-Dx/sqrt(Dx^2+2^2+1)*h*s,1+2/sqrt(Dx^2+2^2+1)*h*s,f(x,1)+h*s/sqrt(Dx^2+2^2+1)], s=0..1, x=-1..1, color=cyan, style=surface):
C4:=plot3d([x-Dx/sqrt(Dx^2+2^2+1)*h*s,-1-2/sqrt(Dx^2+2^2+1)*h*s,f(x,-1)+h*s/sqrt(Dx^2+2^2+1)], s=0..1, x=-1..1, color=cyan, style=surface):
plots:-display(A, B, C1, C2, C3, C4, axes=none);

Exactly the same method can be used for a surface defined parametrically, only the formulas for the normals will be different. For a surface specified implicitly, some technical difficulties arise in the implementation of this method, but I can roughly imagine how they can be solved. Soon I am going to write a post with a procedure to automate all the steps in this plotting.

thick.mw

We get 2 solutions:

restart; 
with(Student:-MultivariateCalculus):
A := [0, 0, 0];
B := [c, 0, 0];
S := [x, 0, z];
solve([Distance(S, A) = a, Distance(S, B) = b], [x, z], explicit) assuming a > 0, b > 0, c > 0;

Obviously, the solutions have a geometric meaning if the triangle inequalities  are hold for the parameters  a, b, c . 

Of course, degenerate cases are also possible when the solution is unique ( z=0  and the points  A, B, S  lie on the same line ). It will be if  a=b+c  or c=a+b  or b=a+c :

Sol:=solve([Distance(S, A) = a, Distance(S, B) = b], [x, z], explicit) assuming a > 0, b > 0, c > 0;
simplify(eval(z,Sol[1]), {a=b+c});
simplify(eval(z,Sol[1]), {b=a+c});
simplify(eval(z,Sol[1]), {c=a+b});

the selection from first level operands (as you wanted):

restart;
expr:=x+sin(x)+ln(y)+10+ln(x+y)^2;
{op(expr)}; 
select(t->type(t,{ln(anything),ln(anything)^anything}), %);

                           {ln(x+y)^2, ln(y)}

Edit.

Maybe you should abandon a polar grid and do as in the example below:

restart;
plots:-polarplot(t,t=0..Pi,color=black, thickness=3, coordinateview = [0 .. Pi, 0 .. Pi]);
plot(t,t=0..Pi,color=black, thickness=3, coords=polar, scaling=constrained, gridlines);

Using the data for the solution  C(x,t) , calculate the derivatives numerically at a point  in which  diff(C(x,t),x,x)<>0  and then  d=diff(C(x,t),t)/diff(C(x,t),x,x)

I replaced  I(t)  with  P(t) , removed 2 strange options 0 and reduced the ranges for  P  and  Q  variables for better viewing:

s0 := 3*10^5:
d := 10^(-3):
delta := 10^4:
b := 5*10^6:
lambda := 4.16:

DEtools:-DEplot([diff(P(t), t) = s0 + P(t)*(-d - delta*Q(t)/(b + Q(t))), diff(Q(t), t) = -lambda*Q(t)], [P(t), Q(t)], t = 0 .. 10, P = 0 .. 20,  Q = 0 .. 20, dirfield = 400, arrows = smalltwo, number = 2, [[0, 4, 0.1], [0, 0.2, 4.1], [0, 7, 0.2], [0, 0.2, 7]], color = red, linecolor = blue, numsteps = 100);

expr:=(x^2 - 1)/x^2*y + (x^2 - 1)/x^2;
subs(x^2 - 1 = Z*x^2, expr);

 
Edit.
 

f := a*x+b;
g := unapply(f, a,b,x);
g(a,b,x); 
g(2,3,x); 
g(2,3,t);

I think the following calculation explains this difference:

restart;
ode1:=diff(y(x),x) = exp(x)*sin(x)+y(x)*cot(x);
my_sol:=[sin(x)*(exp(x)+_C1)];
eval(ode1,y(x)=my_sol);

ode2:=diff(y(x),x) = y(x);
my_sol:=[_C1*exp(x)];
eval( ode2,y(x)=my_sol);

That is if  a<>0  then  [a+b] <> a+[b]

Example:

restart;
X:=LinearAlgebra:-RandomMatrix(50,10, generator=-10..10):
Y:=LinearAlgebra:-RandomMatrix(50,10, generator=-10..10):
B1 := plots:-animate(plots:-pointplot, ['X'[round(kk)], 'Y'[round(kk)], symbolsize = 25, symbol = solidcircle, colorscheme = ["valuesplit"]], kk = 1 .. 50, frames=60);

Just apply the  round  command to all the elements of the matrix.

Example:

A:=<1.1,2.5; -3.53,5.7>;
round~(A);

In my opinion, a conceptually simpler code (especially for a beginner) will just go through a double loop on these lists  X  and  Y :

restart;
X := [3, 6, 4, 13]:   Y := [8, 6, 5, 7]:  n:=nops(X):  k:=0:
for i from 1 to n-1 do
for j from i+1 to n do
k:=k+1;  L[k]:=[[[X[i],Y[i]],[X[j],Y[j]]],sqrt((X[i]-X[j])^2+(Y[i]-Y[j])^2)];
od: od:
L:=convert(L, list);
plot(op~(1,L), style= pointline, symbolsize= 20, symbol= solidcircle);

Addition. The code will be 2 times shorter if we use a nested  seq  instead of loops:

restart;
X := [3, 6, 4, 13]:   Y := [8, 6, 5, 7]:  n:=nops(X):  
L:=[seq(seq([[[X[i],Y[i]],[X[j],Y[j]]],sqrt((X[i]-X[j])^2+(Y[i]-Y[j])^2)], j=i+1..n), i=1..n-1)];
plot(op~(1,L), style= pointline, symbolsize= 20, symbol= solidcircle);

Edit.

Try

if  `and`(seq(a[i] <= b[i] , i=1..min(nops(a),nops(b)) ))  then ...  fi;

# Example:
a:=[1,2,3]:
b:=[1,2,4,5]:
if  `and`(seq(a[i] <= b[i] , i=1..min(nops(a),nops(b)) ))  then true else false  fi;

Try  axis = [tickmarks = [7, subticks = 4], thickness = 2]

Edit. At first I did not notice your second plot.  See a workaround below in a corrected file:

Curves_(1)_new.mw

Example:

Download ellipse.mw