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

It is well known that most of the differential equations cannot be solved symbolically, only numerically. To do this, you must provide an initial condition and use the  numeric  option.

Example:

restart;
Sol:=dsolve({diff(y(x), x) = (6*y(x)^5 - 3*y(x)*x^2 - 20*y(x)^3*x)/(-4*x^3 + 30*y(x)^2*x^2 - 30*y(x)^4 + 7*y(x)^6), y(0)=1}, y(x), numeric);
plots:-odeplot(Sol, [x,y(x)], x=0..1, color=red);

                                            Sol:=proc(x_rkf45) ... end proc
                               

Example:

restart;
plots:-polarplot(cos(3*phi*Pi/180), phi=0..360, angularunit=degrees);

The problem is easily reduced to solving the equation in integers:

restart;
isolve(5*x-7=16*y);

                 {x = 11+16*_Z1, y = 3+5*_Z1}

So  x = 11+16*_Z1  is the set of all the solutions,  _Z1  is an integer.

First you must specify  ma1  as a list:

ma1 := [0., .1941703021, .3203871063, .4089371834, .4712881303, .5145114133, .5435036431, .5617715009, .5718586242, .5756277760, .5744585726]:

subs(ma1[1]=1,ma1);
# Or
subsop(1=1,ma1);

[1, 0.1941703021, 0.3203871063, 0.4089371834, 0.4712881303, 0.5145114133, 0.5435036431, 0.5617715009, 0.5718586242, 0.5756277760, 0.5744585726]
[1, 0.1941703021, 0.3203871063, 0.4089371834, 0.4712881303, 0.5145114133, 0.5435036431, 0.5617715009, 0.5718586242, 0.5756277760, 0.5744585726]

Use the  fieldplot  command from the  plots  package. The grid option controls the number of arrows:

fieldplot(<a, b>, x = -100 .. 100, y = -100 .. 100, arrows = SLIM, grid = [10, 10]);

You can use the combinat:-partition command and then the  select  command to remove duplicate parts:

combinat:-partition(15,8):
select(t->nops(t)=nops({t[]}), %);

[[1, 2, 3, 4, 5], [2, 3, 4, 6], [1, 3, 5, 6], [4, 5, 6], [1, 3, 4, 7], [1, 2, 5, 7], [3, 5, 7], [2, 6, 7], [1, 2, 4, 8], [3, 4, 8], [2, 5, 8], [1, 6, 8], [7, 8]]

restart; 
f:=x->x^2-x+3:
g:=x->3*x-5:
solve(f(x)>=0 or f(x)<0);  # Finding the domain of the function f(x)
h:=f-g:
# Finding the range of the h(x)
m:=minimize(h(x), x=-infinity..infinity);
M:=maximize(h(x), x=-infinity..infinity);

So we have:  (-infinity, infinity)  is the domain of  f(x) ,  [4, infinity)  is the range of  h(x) . Here we use the continuity of the function  h(x)  that takes all intermediate values between  m  and  M .

Use the  seq  command for this.
An example:

restart;
n:=7:
assign(seq(a[i]=i, i=1..n));
P:=piecewise(x<a[1],x, seq(op([x<a[k],(-1)^k*sqrt(0.5^2-(x-(k-0.5))^2)+1]),k=2..n), 1);
plot(P, x=-1..n+2, view=[-1..n+2,-1..2], scaling=constrained, size=[800,400]);