As for the statement of the classical problem see wiki for details.

We assume that:  newly born  n  pairs of rabbits are put in a field; each  pair mates at the age of one month, and at the end of their second month they  produce another pair of rabbits with the probability 3/4 ; and rabbits never die, but continue breeding forever. So we come to the recurrent equation  u(k)=u(k-1)+3/4*u(k-2)  with the initial conditions  u(0)=n, u(1)=n

restart;
u:=unapply(rsolve({u(k)=u(k-1)+3/4*u(k-2),u(0)=n,u(1)=n}, u(k)),k);
expand(u(12));
round(%/n);

Thus, the original number of rabbits will increase approximately 97 times.

restart;
dat := <14.52407334|-162.1310124>:
A := 10^(dat[1]/20.):
phi := dat[2]*Pi/180.:
R:= 0.3036:
f := (Zeta,A,phi) -> cos(Zeta) - R * Zeta * sin(Zeta) - exp(-phi * I)/A:
soln := [solve(f(Zeta,A,phi), Zeta)]:
nos := nops(soln):
# soln := Vector(soln);
soln1:=select(s->Re(s)>0 and Im(s)>0, soln);
L:=sort(soln1,key=(s->Re(s)));
real_part:=Re(L[1]);
imaginary_part:=Im(L[1]);

We can numerically parameterize this curve by polar angle and then do any animations, for example:

restart; 
f1 := sqrt((x+0.500000000)^2+y^2)+sqrt((x-0.500000000)^2+(y-1)^2)+sqrt((x-1.050000000)^2+(y+1)^2)+sqrt((x-0.500000000)^2+(y+0.500000000)^2)+sqrt((x-1)^2+(y-1)^2)-7: 
X := t->fsolve(eval(f1, [x = r*cos(t), y = r*sin(t)]), r = 0 .. 20)*cos(t): 
Y := t->fsolve(eval(f1, [x = r*cos(t), y = r*sin(t)]), r = 0 .. 20)*sin(t): 
plots:-animate(plot, [[X, Y, 0 .. a]], a = 0 .. 2*Pi, frames = 60);

The code below works in Maple 2018.2:
 

restart;
F := proc(t)
local A, B, C;
uses plottools, plots;
A:=line([-2,0], [cos(t)-2, sin(t)], color=blue,thickness=3);
B:=line([cos(t)-2, sin(t)], [t, sin(t)], color=blue,thickness=3);
C:=plot(sin(x), x=0..t, view=[-3..7, -5..5],thickness=3);
display(A,B,C);
end proc:

plots:-animate(F,[theta],theta=0..2*Pi,background=plot([cos(t)-2,sin(t),t=0..2*Pi],thickness=3),
scaling=constrained,axes=none);

You can just skip this option:

restart;
p1 := plot([[0,0], [1,1]]);

If you want the plot of a function  f(x)  without labels, then define the function in operator form, for example:

restart;
plot(x->x^2, -2..2);

restart;
f := a*b+a*c+a+b+c;
s:=select(has, {op(f)}, a) minus {a*b};
# or
s:=select(has,{op(f)} minus {a*b}, a);

Output:       {a, a*c}

restart;
eq2 := ln(2*u^2 + u - 1) = -c - 2*ln(x);
sol:=[solve(eq2,u)];
simplify(sol[1]);
applyop(t->sqrt(simplify(t^2)), [2,2], %);

restart;
f:=x->3+(exp(x^2)+exp(1))/(exp(x^2)-exp(1));
P:=plot(f, -5..5, -3..8, color=blue, discont);
g:=x->f(a)-1/D(f)(a)*(x-a); # Equation of normal to f(x) at the point x=a
x0:=0:
M:=f(0);
d:=sqrt(a^2+(f(a)-g(0))^2):
a:=fsolve(g(0)-M=d, a=1..2);
R:=evalf(d);
y0:=evalf(g(0));
# (x-x0)^2+(y-y0)^2=R^2 - Equation of the circle
plots:-display(P, plottools:-circle([x0,y0],R, color=red), scaling=constrained);

We use a numerical solution ( fsolve  instead of  solve ) because the  solve  command cannot find any solution (returns NULL)

restart;
B[0] := (1/24)*x^4/h^4:
B[1] := -(1/24)*(5*h^4-20*h^3*x+30*h^2*x^2-20*h*x^3+4*x^4)/h^4:
B[2] := (1/24)*(155*h^4-300*h^3*x+210*h^2*x^2-60*h*x^3+6*x^4)/h^4:
B[3] := -(1/24)*(655*h^4-780*h^3*x+330*h^2*x^2-60*h*x^3+4*x^4)/h^4:
B[4] := (1/24)*(625*h^4-500*h^3*x+150*h^2*x^2-20*h*x^3+x^4)/h^4:
piecewise(seq(op([x>n*h and x<=(n+1)*h,B[n]]), n=0..4), 0);

Your surface  f(x,y)  is blocking this ellipse a little, and I raised the ellipse a little:

restart;
with(plots):
f := (x,y) -> x^2 + y^2 - 12*x + 16*y;
display(plot3d(f(x, y), x = -9 .. 9, y = -9 .. 9), pointplot3d([[6, -8, f(6, -8)]], color = red, symbol = solidcircle, symbolsize = 18), view = [-4.2 .. 8.2, -8.2 .. 4.2, -100 .. 100], spacecurve([cos(t), sin(t), 1 - 12*cos(t) + 15*sin(t)+1], t = 0 .. 2*Pi, thickness=2, color=red, orientation = [-15, 68, 5]));

As an easier way you can use  plots:-textplot  command instead, which can easily specify any color.

Example:

restart;
P:=plot(sin(x), x=0..2*Pi, thickness=2, size=[900,400]):
T:=plots:-textplot([[3.5,1,"Plot of "],[4,1,sin(x),color=red]], font=[TIMES,18]):
plots:-display(P,T, scaling=constrained);

When you solve the equation  abs(sin(x)+y^2+y+I*x)  for  y , you are assuming x and y are real numbers. But Maple does not know this and simply solves this equation in the complex domain for the variable  y . For real numbers, it is obvious that   sqrt(u^2+v^2)=0  is equivalent to the system  {u=0,v=0}

A := evalc(abs(sin(x)+y^2+y+I*x)); 
[solve]({op([1, 1, 1], A), op([1, 2, 1], A)})

                                         A := sqrt((sin(x)+y^2+y)^2+x^2)
                                         [ {x = 0, y = -1}, {x = 0, y = 0}]

Example:

restart;
plots:-contourplot(abs(sin(x) + y^2 + y + x*I), x = -6 .. 6, y = -5 .. 5, contours = [seq(1/2^n, n = 1 .. 3), seq(2^n, n = 1 .. 5)], numpoints = 10000, thickness = 2, color = black, filledregions = true, coloring=["LightCyan","Red"]);

In Maple 2018.2, you get the expected result. Of course, formally it will only be a right-sided derivative (as indicated by vv). The code is much shorter if you notice that your function is just  ln(2+5*t)  for  t>=0  and use the differentiate operator  D :

eval(diff(ln(piecewise(t=0,2,2+5*t)),t), t=0);
D(t->ln(2+5*t))(0);

                                               5/2
                                               5/2

The following works as expected:

restart;
B := x^4 - 4*x^3/(2 + p) - 6*(p - 1)*x^2/((3 + p)*(2 + p)) - 4*p*(p - 5)*x/((4 + p)*(3 + p)*(2 + p)) - (p - 1)*(p^2 - 15*p - 4)/((5 + p)*(4 + p)*(3 + p)*(2 + p));
for p from -1 to 5000 do
    A[p] := fsolve(B, x, complex);
end do:
ptlist := convert(A,list):
plots:-complexplot(ptlist, x = -1 .. 1.5, y = -0.5 .. 0.5, style = point);