Important is the way f is defined:

• As a procedure. In that case the variable is established in the definition of f
  

  L := D@@3 + D;
  f := t -> a*t^5 + b*t^4 + sin(c*t);
  L(f)(x);
  

• As an expression. In that case you must explicitly mention the variable.

  restart;
  L := proc(f,t)
     diff(f, t$3) + diff(f,t)
  end proc:
  L( a*x^5 + b*x^4 + sin(c*x), x );
  

What do you mean by "minimize eq1 and eq2"?
Two (different) points (a,b), minimizing eq1 and eq2 respectively?
Or one point (a,b) minimizing eq1*eq2 or eq1+eq2 or...?
Perhaps "eq1" is an expression, depending on a and b which you want to minimaize, subject to the constraint "eq2", an equation containing a and b?

Use decimal points instead of comma's.

g :=(t,s) -> 999.842594+6.793952*10^(-2)*t-9.95290*10^(-3)*t^2 + 
   1.1685*10^(-4)*t^3-1.120083*10^(-6)*t^4 + 
   6.536332*10^(-4)*t^5+8.24493*10^(-1)*s-4.899*10^(-3)*t*s +
   7.6438*10^(-5)*t^2*s-8.2467*10^(-7)*t^3*s +
   5.3875*10^(-4)*t^4*s-5.72466*10^(-3)*s^(3/2) +
   1.227*10^(-4)*t*s^(3/2)-1.6546*10^(-6)*t^2*s^(3/2)+4.8314*10^(-4)*s^2:
g(t,x);

Yes, the four numbers represent [x,y,z,f(x,y,z)]. see ?plot/structure - look for ISOSURFACE(A)

Now for removing the last number in each list:

ptslist[1..-1,1..-1,1..-1,1..3]:

there must have been more backquotes than the two you removed

restart;
a `:=` 4;
Error, missing operator or `;`

restart;
`a `:=` 4`;
                             a  :=  4
a, `a `;
                             a,  4
`a := 4`;
                             a := 4
b := %;
                            b :=  a := 4

All your equations contain

diff(phi,x)   or   diff(phi,y)

I suppose that phi is assigned earlier, as an expression containing x and y (otherwise diff(phi,x) simplifies instantly to zero).
Or do you mean

diff(phi(x,y),x)

Please clarify.

If so (i.e: phi is a not specified function of x and y), you can do

restart;
alias(phi=phi(x,y)):
eqns := {-2/phi^3*u[0]^2*diff(phi,x)-4/phi^3*w[0]*diff(phi,x)-2*u[0]/phi^3*diff(phi,x)^2, 
     -6/phi^4*u[0]*diff(phi,x)*v[0]+12*v[0]/phi^4*diff(phi,x)^2, 
     -2/phi^3*v[0]*diff(phi,x)+2/phi^3*w[0]*diff(phi,y)}:
U := solve( eqns, {u[0],v[0],w[0]} );  # a sequence of solution sets
nops( [U] );  # make a list, and count the number of its elements.
                               3

so there are three solutions

Does this what you want?

restart;
f := ln( (r1(x,y)+r2(x,y)+d1)/(r1(x,y)+r2(x,y)-d1) ):
df := simplify(diff(f,x));
r1 := (x,y) -> sqrt((x-x1)^2+(y-y1)^2+z^2): alias(R1=sqrt((x-x1)^2+(y-y1)^2+z^2)):
r2 := (x,y) -> sqrt((x-x2)^2+(y-y2)^2+z^2): alias(R2=sqrt((x-x2)^2+(y-y2)^2+z^2)):
d1 := sqrt((x2-x1)^2+(y2-y1)^2): alias(D1=sqrt((x2-x1)^2+(y2-y1)^2)):
df;

So you really want symbolic solutions?
OK here it comes.
I make all parameters symbolic, there is no need to give some of these as floats.

Hrs := -n*Np*(Ans*r^(n*Np-1)-Bns*r^(-n*Np-1))*cos(n*Np*theta):
mu1 := k1*k2*Hrs^(k2-1):
eq1 := Bno*Ro^(-n*Np-1) = Ans*Ro^(n*Np-1)+Bns*Ro^(-n*Np-1):
eq2 := -Bno*Ro^(-n*Np-1) = mu1*(Ans*Ro^(n*Np-1)-Bns*Ro^(-n*Np-1)):
r := Rs:
Hrs :=  -n*Np*(Ans*r^(n*Np-1)-Bns*r^(-n*Np-1))*cos(n*Np*theta):
mu1 := k1*k2*Hrs^(k2-1):
eq3 := Ans*Rs^(n*Np-1)+Bns*Rs^(-n*Np-1) = Ana*Rs^(n*Np-1)+Bna*Rs^(-n*Np-1):
eq4 := mu1*(Ans*Rs^(n*Np-1)-Bns*Rs^(-n*Np-1)) = Ana*Rs^(n*Np-1)-Bna*Rs^(-n*Np-1):
eq5 := Ana*Rm^(n*Np-1) + Bna*Rm^(-n*Np-1) = 
       Anm*Rm^(n*Np-1) + Bnm*Rm^(-n*Np-1)+M/(mu2*(1-(n*Np)^2)):
eq6 := Ana*Rm^(n*Np-1) - Bna*Rm^(-n*Np-1) = 
  mu2*(Anm*Rm^(n*Np-1) - Bnm*Rm^(-n*Np-1)+M/(mu2*n*Np*(1-(n*Np)^2))) - 
   Mrad/(n*Np):
eq7 := Anm*Rr^(n*Np-1) + Bnm*Rr^(-n*Np-1) + M/(mu2*(1-(n*Np)^2)) = Anr*Rr^(n*Np-1):
eq8 := mu2*(Anm*Rr^(n*Np-1)-Bnm*Rr^(-n*Np-1)+M/(mu2*n*Np*(1-(n*Np)^2))) - Mrad/(n*Np) = mu3*Anr*Rr^(n*Np-1):

We see that only eq2 and eq4 are nonlinear, so we make a linear system for the other six equations:

with(LinearAlgebra):
A,b := GenerateMatrix( {eq1, eq3, eq5, eq6, eq7, eq8}, 
  [Ana, Anm, Anr, Ans, Bna, Bnm, Bno, Bns] ):

A is 6 ×8-matrix. When we solve, we get a solution vector with two free parameters:

sols := LinearSolve(A,b,free='alpha'):

The indexing of the α_i's is unpredictable, so we recover

awa := [seq(has~(sols,alpha[i]),i=1..8)]:
aind := NULL: 
for i to 8 do 
  if convert(convert(awa[i],list),`or`) then 
  aind := aind,i 
  end if 
end do:

Now substitute in the two nonlinear equations, and solve for the free variables.

s1 := Equate( Vector([Ana, Anm, Anr, Ans, Bna, Bnm, Bno, Bns]), sols ):
salpha := solve ( subs( s1, {eq2,eq4} ), {seq(alpha[i],i=aind)} );

          [Length of output exceeds limit of 1000000]

Voila

A more elegant way is to make the function continuous in the origin by putting the value for (0,0) in a remember table:

f := (x,y) -> sin(sqrt(x^2+y^2))/sqrt(x^2+y^2): f(0.,0.) := 1:

Caution: not f(0,0) becasuse grid points are floats, so plot3d is looking for an index [0.,0.] (floats), and not for [0,0] (integers).

plot3d( f, -3..3, -3..3 );

or

plot3d( 'f(x,y)', x=-3..3, y=-3..3 );

result in a graph with no hole.
Caution:

plot3d( f(x,y), x=-3..3, y=-3..3 );

doesn't work, because f(x,y) is immediately replaced by the expression sin(sqrt(x^2+y^2))/sqrt(x^2+y^2) in the plot command (premature evaluation!).

Maple makes a plot bij calculating f(x,y) in a number of gridpoints, derived from the ranges for x and y.
In this case f(0,0) is not defined ("division by zero"), so we have to prevent maple from taking (0.0) as a grid point. An easy way to do this, is to shift the ranges by a small amount

plot3d( sin(sqrt(x^2+y^2))/sqrt(x^2+y^2), x=-3.01..2.99, y=-3.01..2.99);

g := x -> piecewise( x<3, 2*x, 2*x+5 ):
plot(g(x));

Not exactly the same, but perhaps ?member or ?subset will do:

a := {1,2,3}:
A := {1,2,3,4}:
member( 1, A );
                              true
a subset A;
                              true

By the assignment

A[m] := ((Ck[m](i,j))/((lam(i,i)-lam(j,j))-(2*m));

you make A[m] a number, not a matrix
If A is a table of zero-matrices, they can be filled by

 A[m](i,j) := some number 

Additionally: count parentheses (and remove some)!

Perhaps yoy can use Markiyan Hirnyk's answer to this previous question.

I'm afraid that there is no other possibiliy than to tell Maple that the two variables are interdependent (btw, what does dx dy mean in this case?) You could do that by

c3 := Int(Int( subs( [x=f(x,y), y=g(x,y)], x*y+y ), x), y);