In fact, the  plot is shown incorrectly  and  pp(2020) is calculated incorrectly . The coefficients of the polynomial differ too much from each other. Simply increasing  Digits  will not help here. It is unclear to what extent  Digits  should be increased. Conversion of all constants into exact fractions radically solves the problem.. Then all calculations are done exactly:

restart;
a := 0.340699639252428*10^13;
b := -0.118229737138742*10^11;
c := 0.175816922773262*10^8;
d := -14523.7138112711;
e := 7.19782673654200;
f := -0.00214008825847444;
g := 0.353463615941623*10^(-6);
h := -0.250170509163729*10^(-10);
a,b,c,d,e,f,g,h:=convert~([a,b,c,d,e,f,g,h],fraction)[];

pp := x -> a + b*x + c*x^2 + d*x^3 + e*x^4 + f*x^5 + g*x^6 + h*x^7;

plot(pp(x), x = 1960 .. 2100, gridlines = true, size = [800, 400]);
pp(2020);
evalf(%);

Now all the contradictions between the values of the function  pp  on the plot and the values of the function  pp  calculated at separate points  ( like pp(2020)) disappear.

Example:

restart;
myTable:=table(["a"=1,"b"=-1,"c"=1]);
select(i->myTable[i[]]=1,{indices(myTable)});

                                               {["a"], ["c"]}

Edit.

If we define  f  as a procedure, then we get a more convenient and compact syntax for specific calculations: 

restart;
with(orthopoly):
f:=unapply(convert(diff( P(n, cos(theta)), theta, theta),diff), n):
int( f(2), theta=0..Pi/4);
int( f(3), theta=0..Pi/4);

As an alternative you can use  diff  instead of  D :

restart;
g:=x->x^(1/3);
f:=x->convert(g(x),surd);
diff(f(x),x);

It's very simple:

restart;
A:=[1,2,6,7,9,10,15,17]:
A[3..5];

                                             [6, 7, 9]

restart;
f:=(a*x1+b*x2)/x3+c*x1^2/(e*x2+d);
normal(f);
numer(%);

Or

restart;
(a*x1+b*x2)/x3+c*x1^2/(e*x2+d)=f;
normal(%);
%*denom(lhs(%));

Code in Maple 2018.2 (I have no newer versions):

restart;
B:=sqrt( (-4*u^(1/3)+1)*u^(4/3));
A:=1/(-12*u+3*u^(2/3)-3*B);
res:=Int(A,u);
res1:=IntegrationTools:-Change(res,t=u^(1/3));
value(res1);
simplify(%) assuming t>0;
eval(%, t=u^(1/3)); 

You can get several more solutions from one solution if you use transformations that transform a cubic integer lattice into itself. They will be, for example, reflections about the coordinate planes, rotations around the coordinate axes at angles that are multiples of 90 degrees, and some others. But in terms of geometry, you end up with the same pair of vectors. New solutions can be easily obtained using the brute force method. First we find the list L and then from the list  L, the list  L1  is obtained, where none of the coordinates is equal to 0:

restart;
k:=0:
for a from 0 to 9 do
for b from a to 9 do
for c from b to 9 do
for x from -9 to 9 do
for y from -9 to 9 do
for z from -9 to 9 do
if 4*(a*x+b*y+c*z)^2=3*(a^2+b^2+c^2)*(x^2+y^2+z^2) and a^2+b^2+c^2<>0 and x^2+y^2+z^2<>0 and a*x+b*y+c*z>0 then k:=k+1;
L[k]:=[<a,b,c>,<x,y,z>] fi;
od: od: od: od: od: od:
L:=convert(L,list):  
nops(L);
L1:=select(t->`and`(seq(t[1][i]<>0,i=1..3),seq(t[2][i]<>0,i=1..3)),L):
nops(L1);

The very first pair of vectors in the lists  L  and  L1  is your original example, in which the first and third coordinates are reversed (reflection of space relative to the plane  x=z ). To reduce the search, I chose the vector  <a,b,c>  with non-decreasing coordinates.

This can be done in many ways. Here's one:

A := [1,0,3,4,5,5]: B := [0,1,2,6,3,6]:
nops(select(`>`, A - B, 0));

                                           3

The simple recursive procedure does the work:

restart;
a:=proc(n)
local m;
option remember;
if n=1 then return 1 elif n=2 then return 0 else
for m from 1 to n-2 do
if a(n-1)=a(n-1-m) then return m fi; od;
min(seq(abs(a(n-1)-a(k)),k=1..n-2)) fi;
end proc:

Examples of use:

seq(a(n), n=1..100);
plots:-listplot([%], style=point, symbol=solidcircle, color=red);

Edit. The procedure above works about 2 times faster than tomleslie's method.

In real domain the two-argument function  arctan(y,x)  returns the polar angle  alpha  of any point with the coordinates  (x,y)  on the plane in the range  -Pi < alpha <= Pi . Below are all the variants as it expressed through the standard function arctan(x) :

arctan(y,x) assuming x>=0,y>=0;
arctan(y,x) assuming x>=0,y<0;
arctan(y,x) assuming x<0,y>=0;
arctan(y,x) assuming x<0,y<0;

If I understand your problem correctly, then you want to get derivatives of different orders from the original expression  u[0]  using the formulas for differentiating the sum and the product at each step of the for-loop. See the following code for your specific example:

restart;
u[0] := exp(-x^2);
n:=10:
d := r->diff(coeff(r,u[0]),x)*u[0]+coeff(r,u[0])*diff(u[0],x):

for k from 1 to n do 
if type(u[k-1],`+`) then u[k] := map(d,u[k-1]) else 
u[k] := d(u[k-1]) end if; 
print(u[k]); 
end do:

Sx:=1/sqrt(2)*Matrix([[0,1,0],[1,0,1],[0,1,0]]);
lam,v:=LinearAlgebra:-Eigenvectors(Sx);
map(r->Vector(r/~LinearAlgebra:-VectorNorm(Vector(r),2)),convert(v^%T,listlist)); 

The problem is easy to solve if you just cut out regions with asymptotes from the plot:

restart;
func:=unapply(-tan((-4*t+x+2*y)/(2*sqrt(15)))/(2*sqrt(2)), x,y,t):
y1:=solve((-4*t+x+2*y)/(2*sqrt(15))=-Pi/2,y);
y2:=solve((-4*t+x+2*y)/(2*sqrt(15))=Pi/2,y);
y3:=solve((-4*t+x+2*y)/(2*sqrt(15))=3*Pi/2,y);
y4:=solve((-4*t+x+2*y)/(2*sqrt(15))=5*Pi/2,y);
A:=plot3d(func(x,y,4), x=-10..10, y=eval(y1+0.1..y2-0.1,t=4), style=surface):
B:=plot3d(func(x,y,4), x=-10..10, y=eval(y2+0.1..y3-0.1,t=4), style=surface):
C:=plot3d(func(x,y,4), x=-10..10, y=eval(y3+0.1..y4-0.1,t=4), style=surface):
plots:-display(A,B,C, scaling=constrained);

restart;
eq := cos(2*x)-sin(x)*(sin(x)^2+1)^(1/2)+cos(x)^2*sin(x)/(sin(x)^2+1)^(1/2):
solve(subs(x=arcsin(y), eq));