x0,y0,z0:=1,2,3:
eq:= u*~[cos(t),sin(t),0] + (1-u)*~[x0,y0,z0]:
plot3d( eq, u=0..1,t=0..2*Pi, scaling=constrained );

For the 2nd order derivative you must use

fdiff(dydata, [1,1], 0);

[Note that to compute this it would be much better to use the RHS of the ODE].

With your new constraints Maple cannot find a feasible integer point because such point does not exist. i.e. your constraints are incompatible.
You can see this using Maple (but also by hand):

The constraints are:

m := 9: n := 9: p := 3:
constr:= seq(add(Y[i,j],i=1..m)=1,j=1..n),
         add(X[i],i=1..m)=p,
         seq(seq( Y[i,j]<=X[j],i=1..m),j=1..n);

By symmetry we may take
X[1]:=1;X[2]:=1;X[3]:=1; for i from p+1 to m do X[i]:=0 od;

But Y[i,j]<=X[j] implies:

for i to m do for j from p+1 to n do Y[i,j]:=0  od: od;

and evaluating
constr;  # ==>

Y[1, 1]+Y[2, 1]+Y[3, 1]+Y[4, 1]+Y[5, 1]+Y[6, 1]+Y[7, 1]+Y[8, 1]+Y[9, 1] = 1, Y[1, 2]+Y[2, 2]+Y[3, 2]+Y[4, 2]+Y[5, 2]+Y[6, 2]+Y[7, 2]+Y[8, 2]+Y[9, 2] = 1, Y[1, 3]+Y[2, 3]+Y[3, 3]+Y[4, 3]+Y[5, 3]+Y[6, 3]+Y[7, 3]+Y[8, 3]+Y[9, 3] = 1, 0 = 1, 0 = 1, 0 = 1, 0 = 1, 0 = 1, 0 = 1, 3 = 3, Y[1, 1] <= 1, Y[2, 1] <= 1, Y[3, 1] <= 1, Y[4, 1] <= 1, Y[5, 1] <= 1, Y[6, 1] <= 1, Y[7, 1] <= 1, Y[8, 1] <= 1, Y[9, 1] <= 1, Y[1, 2] <= 1, Y[2, 2] <= 1, Y[3, 2] <= 1, Y[4, 2] <= 1, Y[5, 2] <= 1, Y[6, 2] <= 1, Y[7, 2] <= 1, Y[8, 2] <= 1, Y[9, 2] <= 1, Y[1, 3] <= 1, Y[2, 3] <= 1, Y[3, 3] <= 1, Y[4, 3] <= 1, Y[5, 3] <= 1, Y[6, 3] <= 1, Y[7, 3] <= 1, Y[8, 3] <= 1, Y[9, 3] <= 1, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0

As you can see, there are some "0=1" here !

You will have to reconsider the constraints.
 

Only a few corrections were needed.

p-Median-ok.mw

This should work.

plots[fieldplot3d]([2*x, 2*y, 2*z], x = -1 .. 1, y = -1 .. 1, z = -1 .. 1, grid = [5, 5, 5], arrows = SLIM,
colorscheme=["xyzcoloring", (x, y, z)->sin(x)*sin(y)*sin(z)]);

It is accepted syntactically but it seems to be ignored.
In the help:
"Most but not all plotting commands support this option".

They are not quite similar. The multi-ellipses are always convex while the lemniscates could be even nonconnected.
On the other side, when the radius tends to infinity, both become almost circles.

lemni:=proc(p,R:=1,M:=10)
local x,y, z:=indets(p,symbol)[1], P:=eval(simplify(p),z=x+I*y):
P:=evalc(P*conjugate(P)-R):
plots[implicitplot](P, x=-M..M,y=-M..M, numpoints=10000)
end:

lemni((z^7-1)/(z-1),1);
lemni((z^7-1)/(z-1),1000);

M:=LinearAlgebra:-JordanForm(A, output='Q');

(F must be even.)

alpha := (1/8)*Pi:  R := 100: Wi := 3/10: H := 0: 
Eq1 := diff(F(eta), eta, eta, eta)+2*alpha*R*F(eta)*(diff(F(eta), eta))+4*alpha^2*(diff(F(eta), eta))+Wi*(8*alpha*F(eta)*(diff(F(eta), eta, eta, eta))+32*alpha^3*F(eta)*(diff(F(eta), eta)))-H*R*(cos(alpha*eta)^2*(diff(F(eta), eta))-alpha*F(eta)*sin(2*alpha*eta)) = 0:
IC1 := F(0) = 1, F(-1) = 0, F(1) = 0:
MyIC1:= F(0) = 1, D(F)(0) = 0, F(1) = 0:
MyIC2:= F(0) = 1, D(F)(0) = 0, F(-1) = 0:
smy1:=dsolve({Eq1, MyIC1}, F(eta), numeric):
smy2:=dsolve({Eq1, MyIC2}, F(eta), numeric):
p1:=plots:-odeplot(smy1):
p2:=plots:-odeplot(smy2):
plots:-display(p1,p2);

Explore could be helpful. It is not clear if a solution exists.
 

s:=dsolve({Eq1, F(0)=1, D(F)(0)=p, (D@@2)(F)(0)=q}, F(eta), numeric, parameters=[p,q]);
P:=proc(p,q)
  s(parameters=[p,q]);
  plots:-odeplot(s, -1..1)
end;
Explore( P(p,q), p=-10. .. 10., q=-10. .. 10., initialvalues=[p=0,q=0] ); 

Edit: forgot the interval in odeplot

Set
Digits:=45;
and it works. (40 is not enough.)
t0 := 0.551325384496649794326454673778951177891621325e-1

The question is why do you need Analytic for a simple equation in x
( cos(Pi*x*p)=q).

`index/defseq1` := proc(L)   local n:=op(L);  
      n^2   end proc:
a:=table(defseq1);

a[10] + a[n+1];
   100+(n+1)^2

You are not in a Boolean context, so you must use:

(t->evalb(t="."))(9);

You have used FDM in another interval.
Setting N=5 the two methods will agree.
If you keep N=10, you must set h:=1 but then fsolve cannot solve the system;  try to use DirectSearch instead.

If you want the big picture of the solutions of your system, of even if you want all the solutions in a prefered order, use Groebner bases.

For example, a system with 3 unknowns

P1(x,y,z,a)=0, P2(x,y,z,a)=0, P3(x,y,z,a)=0.

can usually be rewritten in an equivalent form:

Q1(x,a)=0, Q2(x,y,a)=0, Q3(x,y,z,a)=0.
and here the structure is clear (Pi,Qi are polynomials).
See ?Groebner
But as Carl said, we cannot treat an "abstract" case without an example, it would imply a course in Polynomial algebras, ideals etc.