When simplifying complex expressions, it is often useful to replace some subexpression with a name and use subs  or  algsubs  commands, and instead of  simplify  use the  factor  command with the required options:

The set of points  z=x+I*y  on the complex plane satisfying this equation  Eq  is empty . See below

restart;
z:=x+I*y:
Eq:=evalc(abs(2*z-I)=Im(z+1-I)); # The equation Eq is equivalent to the system S
S:={Eq^2, y>=1};
solve(Eq);
is(-3*y^2+2*y<0) assuming y>=1;

5 points define a polynomial of degree <=4 :

CurveFitting:-PolynomialInterpolation([[1,2],[2,3],[3,4],[4,6],[5,10]], x);
plot(%, x=0..5, tickmarks=[6,11], gridlines, scaling=constrained, size=[300,500]);

Addition. If you still want to get a polynomial of degree 5, then you need to specify one more point. I have defined it with a parameter  a . We see that as a parameter we can take any  a<>17 :

CurveFitting:-PolynomialInterpolation([[1,2],[2,3],[3,4],[4,6],[5,10],[6,a]], x);

In the solution above, nowhere are the forces 70 lb and 60 lb taken into account. 
In this problem, the unknown quantities are: the reaction of the support  OA (say F), the weight of the pot (say P) , and the 1st coordinate  x  of point  A. To find them, we compose a system of 3 equations with 3 unknowns. After finding  x , the coordinates of all points will be known and it will be easy to find all the angles. The sum of all forces applied to point  A  should equal 0.

restart;
local O;
O:=<0,0,0>: A:=<x,6,2>: B:=<-3,0,4>: C:=<5/2,0,5*sqrt(3)/2>:
AC:=C-A; AB:=B-A; AO:=O-A;
assume(x>0): 
vOA:=A/sqrt(A.A)*F; # Force vector in the direction of OA 
vAB:=AB/sqrt(AB.AB)*70; # Force vector in the direction of AB   
vAC:=AC/sqrt(AC.AC)*60; # Force vector in the direction of AC  
vP:=<0,0,-P>; # Force vector vertically down 
Sys:={seq((vOA+vAB+vAC+vP)[i]=0, i=1..3)};
Sol:=fsolve(Sys, {F,P,x});
x:=eval(x,Sol);
`Angle BAC`=evalf(arccos(AB.AC/sqrt(AB.AB)/sqrt(AC.AC))*180/Pi)*degrees; 
`Angle OAC`=evalf(arccos(AO.AC/sqrt(AO.AO)/sqrt(AC.AC))*180/Pi)*degrees;

Edit.

i:=1: k:=3:
G[i,k]=abs(u[i] - u[k])-(u[i]-u[k])^2;

You can generate your table as a matrix:

Matrix(5, (i,j)->`if`(i>=j,'y[i,j]'=y[i,j], ``));

restart;                    
T := (p*a^(-Phi(xi))+q+r*a^Phi(xi))/ln(a):                    
u[0] := C[0]+C[1]*a^Phi(xi)+C[2]*a^(2*Phi(xi)):                          
u[1] := diff(u[0], xi):               
d[1] := C[1]*a^Phi(xi)*T*ln(a)+2*C[2]*a^(2*Phi(xi))*T*ln(a):         
u[2] := diff(d[1], xi):     
d[2] := C[1]*a^Phi(xi)*T*ln(a)*(p*a^(-Phi(xi))+q+r*a^Phi(xi))+C[1]*a^Phi(xi)*(-p*a^(-Phi(xi))*T*ln(a)+r*a^Phi(xi)*T*ln(a))+4*C[2]*a^(2*Phi(xi))*T*ln(a)*(p*a^(-Phi(xi))+q+r*a^Phi(xi))+2*C[2]*a^(2*Phi(xi))*(-p*a^(-Phi(xi))*T*ln(a)+r*a^Phi(xi)*T*ln(a)):                                     
expand((2*k*k)*w*beta*d[2]-(2*alpha*k*k)*d[1]-2*w*u[0]+k*u[0]*u[0]):     
subsindets(%,identical(a^Phi(xi)), t->z);
collect(%, z);
subs(z=a^Phi(xi), %);

Use  unapply  instead of a function in arrow notation:

restart;
K := x^3*y^4 + 6*x^2*y^3 + 3*x*y^4 + x^2*y^2 + 2*x*y^3;
f := unapply(K, x,y); 
f(t, x);

You can see for yourself that when using arrow notation  f:=(x,y)->K , the right-hand side  K  must be explicitly specified in terms of  x  and  y .

restart;
V:= Vector([3 * x[1] ^ 4 * x[2] ^ 2 + 7 * x[1] ^ 3 * x[2] ^ 3 + 6 * x[1] ^ 2 * x[2] ^ 4 + 7 * xi[1] * xi[2] ^ 5 + 65 * xi[1] ^ 3 + 76 * xi[1] ^ 2 * xi[ 2] + 56 * xi[1] ^ 2 * xi[3], 13 * x[1] ^ 6 - 7 * x[1] ^ 5 * x[2] + 30 * x[1] ^ 4 * x[2] ^ 2 + 75 * x[1] ^ 3 * x[2] ^ 3 + 130 * x[1] ^ 2 * x[2] ^ 4 + 54 * xi[1] * xi[2] ^ 5 + 43 * xi[2] ^ 6 + xi[2] ^ 2 * xi[3] + 76 * xi[2] ^ 2 * xi[4] + 43 * xi[2] * xi[3] ^ 2 ]);
 
map(v->select(t->degree(t)=3, v), V);
       

To get the list of coefficients use the  coeffs  command.

You can use the  Student:-LinearAlgebra:-RotationMatrix  command to parameterize a circle in 3D. In the example below  M  is rotation matrix around the axis defined by the vector  <2,3,-4>  by the angle  t ,  A  is the center of the circle,  B  is the point on this circle: 

restart;
M:=Student:-LinearAlgebra:-RotationMatrix(t, <2,3,-4>);
A:=<1,2,3>: 
B:=<2,4,5>:
Eq:=A + M.(B-A); # The parametric equations of the circle
plots:-spacecurve(Eq, t=0..2*Pi, color=red, axes=normal, scaling=constrained, view=[-2..5.5,-2..5.5,-2..5.5]);

 
Note that the vectors  <2, 3, -4>  and  B - A  are perpendicular.

restart;
y := t -> 0.7*sin(300*t);
filenm := "https://www.gw-openscience.org/GW150914data/P150914/fig2-unfiltered-template-reconstruction-H.txt";
K := ImportMatrix(filenm, source = Matlab, skiplines = 1);
R:=[seq([K[i,1],abs(K[i,3]-y(K[i,1]))], i=1..3441)]:
plot([K[() .. (), [1, 3]], y(t)], t = 0.35 .. 0.4, color=[red,blue]);
plot(R, color=green, view=[0.35..0.4,0..1]);

I understood the plot of residuals as the plot of the modulus of the difference of two plots (blue and red).

Edit. If you just take the differences (no modulus), you get a graph like Carl's.

R:=[seq([K[i,1],y(K[i,1])-K[i,3]], i=1..3441)]:
plot(R, color=green, view=[0.35..0.4,-1.5..1.5]);

When you use  break, you only get out of one innermost loop (the rest continue to work).
To solve the problem, just replace two lines of code 
print(dedges);
break;
with one line
return dedges;

Edit. For any number of pairs of edges, you can always limit yourself to one cycle. Carl has already answered your second question. Below is another method that does not use the Iterator package. This method is suitable also for older versions of Maple:

restart;
with(GraphTheory):
with(SpecialGraphs):
g:=CompleteGraph(6);
g1:=DeleteEdge(g, {3,2}, inplace = false);

Conremovedges:= proc(G::Graph,Paircrossedges)
local n, P, p, dedges;
uses combinat;
n:=nops(Paircrossedges);
P:=permute([1$n,2$n],n);
for p in P do
dedges:={seq(Paircrossedges[i][p[i]], i=1..n)}; 
if EdgeConnectivity(DeleteEdge(G, dedges, inplace = false))=3 then
return dedges;
end if;
end do:
"not found";
end proc:

# Example of use:
Paircrossedges:=[[{1,6},{3,5}],[{2,5},{1,4}],[{2,6},{3,4}]];
Conremovedges(g1,Paircrossedges);

If you want the square root not to be grayed out when using the inert form, but to look as usual, then use the  InertForm  package with the option  inert=false :

InertForm:-Display(1/%sqrt(2), inert=false);

If you just need a procedure that returns the distance from a specific vertex  v  to a specific edge  e  in a graph  G , then here it is:

restart;
Dist:=proc(v,e,G)
uses GraphTheory;
min(Distance(G,v,e[1]), Distance(G,v,e[2]));
end proc:

Examples of use:

G := GraphTheory:-SpecialGraphs:-PetersenGraph():
E := convert(GraphTheory:-Edges(G),list);
Dist(1,{8,9},G); # Distance from the vertix 1 to the edge {8,9}
Dist~(1,E,G); # List of distances from vertex 1 to all edges from the list E 

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.