#
# Extract the parameter values from the DirectSearch and
# LSSolve results tables. Put these into matrices, together
# with an "index" row and an "index" column.
#
  LSMat:= Matrix( max([indices(ansLS, 'nolist')])+1,
                  numelems(ansLS[1][2])+1
                ):
  LSMat[1,1]:= "n":
  LSMat[1,2..-1]:= convert~( < lhs~(ansLS[1][2])[] >,
                             string
                           ):
  LSMat[2..-1,1]:= < seq
                     ( j,
                       j=1..n-1
                     )
                   >:
  LSMat[2..-1, 2..-1]:= < seq
                          ( Transpose
                            ( < rhs~(ansLS[j][2]) > ),
                            j=1..n-1
                          )
                        >:
  DSMat:= Matrix
          ( max
            ( [ indices
                ( ansDS, 'nolist')
              ]
            )+1,
            numelems
            ( ansDS[1][3] )+1
          ):
  DSMat[1,1]:= "n":
  DSMat[1,2..-1]:= convert~( < lhs~(ansDS[1][3])[] >,
                             string
                           ):
  DSMat[2..-1,1]:= < seq
                     ( j,
                       j=1..n-1
                     )
                   >:
  DSMat[2..-1, 2..-1]:= < seq
                          ( Transpose
                            ( < rhs~(ansDS[j][3]) > ),
                            j=1..n-1
                          )
                        >:
#
# So are the results from the LinearSolve() and
# DirectSearch() processes very different?? They
# shouldn't be! For amusement, find the biggest
# absolute difference for any parameter
#
   mi:= max[index](DSMat-LSMat);
  (DSMat-LSMat)[mi];
#
# Differences probably too small to be noticeable
# when data is "rounded" for export to Excel
#
# Send data to two different Excel spreadsheets,
# one for the LSSolve() results and one for the
# DirectSearch() results.
#
# OP will have to change 'fname' to something
# appropriate for his/her machine, OS
#
  fname:= "C:/Users/TomLeslie/Desktop/":
  ExcelTools:-Export( LSMat, cat(fname, "OPTLSres.xlsx"), 1, "B2"):
  ExcelTools:-Export( DSMat, cat(fname, "OPTDSres.xlsx"), 1, "B2"):

#############################
# Initialize and get version
#
  restart;
  kernelopts(version);
#
# OP's code
#
  imp_fun := -4*x + 10*(x^2)*(y^(-2)) + y^2 =11;
  c := 2;
  s:= evalf( solve( subs( x = c, imp_fun)));
  m1 := evalf( subs ( { x = c, y = s[1] }, implicitdiff( imp_fun, y,x)));

#############################
# How I'd (probably?) do it.
#
# Initialize and get version
#
  restart;
  with(plots):
  kernelopts(version);
  xVal:=2:
#
  imp_fun := x^3+y^3 = 9*x*y;;
  s:= [ evalf
        ( solve
          ( eval
            ( imp_fun,
              x=xVal
            )
          )
        )
      ];
  m:= [ seq
        ( evalf
          ( eval
            ( implicitdiff
              ( imp_fun, y, x
              ),
              [ x=xVal, y=s[j] ]
            )
          ),
          j=1..numelems(s)
        )
      ];
#
# Plot the curve and the solution points for x=2
#
  cols:=[ red, green, blue, black]:
  display
  ( [ implicitplot
      ( imp_fun,
        x=-5..5,
        y=-5..5,
        color=black,
        gridrefine=4
      ),
      pointplot
      ( [ seq
          ( [xVal, s[j]],
            j =1..numelems(s)
          )
        ],
        color=[ red, green, blue],
        style=point,
        symbol=solidcircle,
        symbolsize=20
      ),
      plot
      ( [ seq
          ( s[j]+m[j]*(x-xVal),
            j=1..3
          )
        ],
        x=-5..5,
        color=[ red, green, blue]
      )
    ]
  );

  restart;
  with(plots):
#
# Define the data: life would be so much easier if the
# OP could be bothered to do this. It would mean that I
# don't have to spend time generating "artificial" data
#
  g:=2926: nL:=26:
  X:=<seq(j, j=0..evalf(2*Pi),evalf(2*Pi/g))>:
  M:=Matrix( nL, g, (i,j)->i/100*sin(X[j])):
#
# Now that data has been generated, which the OP should
# have provided in the first place, then do a simple
# plot of all the data
#
  display( [seq(pointplot( X, M[j,..], style=line), j=1..nL)]);

  tol:= 1e-06;
#
# Define value for maximum number of iterations
#
  n:= 10;
#
# Inuitialze array for results of successive iterations
#
  x:= Array(0..n):
  x[0]:= 0.1;
#
# Define a couple of functions
#
  f:= z -> z^2-3;
  h:= z -> 2*z;
#
# Perform successive iterations
#
  for k from 1 to n do
      x[k]:= evalf(x[k-1]-f(x[k-1])/h(x[k-1]));
      if   abs(x[k]-x[k-1]) < tol
      then print("Number of iterations"=k);
           print("approximate solution"=x[k]);
           print( f(x[k]) );
           break;
      end if;
   # x[k-1]:= x[k]; this line does nothing useful so comment out
  end do:
#
# Display values of successive iterations
# (allowing for Maple's reluctance to display
# large amounts of data whan this isn't necessary)
#
  interface(rtablesize=n+1);
  x();
  interface(rtablesize=10);
#
# Produce plot
#
  plot( [f(p), f(x[k-1])+(p-x[k-1])*h(x[k-1])], p=-1..3, colour=[red,black]);
 

  restart;
  kernelopts(version);
  with(plots):
  fcns := {T(eta), f(eta)}:
  ep:= .1: M := 1: kp := .5: n := 1:
  ec:= .1: pr := 1: s := .1: N := 5:
  sys:= diff(f(eta),eta$3)+f(eta)*diff(f(eta),eta$2)-diff(f(eta),eta)*diff(f(eta),eta)+ep*ep+(M+1/kp)*(ep-diff(f(eta),eta)) = 0,
        diff(T(eta),eta$2)+pr*(f(eta)*diff(T(eta),eta)-n*diff(f(eta),eta)*T(eta))+pr*(ec*diff(f(eta),eta$2)*diff(f(eta), eta$2)+ec*(M+1/kp)*diff(f(eta),eta)^2+s*T(eta)) = 0:
  bc:= f(0)=0, D(f)(0)=1, D(f)(N) = ep, T(0) = 1, T(N) = 0:
  R:= dsolve(eval({bc, sys}), numeric, output=listprocedure):
  psi := [-0.9e-1, -0.7e-1, -0.4e-1, -0.1e-1, 0, 0.1e-1, 0.4e-1, 0.7e-1, 0.9e-1]:
  for i to 9 do
      plt[i]:=odeplot( R, [eta, psi[i]/f(eta)], eta=0..5, axes=boxed);
  end do:
  display( [seq(plt[i], i=1..9)]);

  cols:= [ "Purple", "Maroon", "Red", "Orange",
           "Gold", "Yellow", "Green", "LightSkyBlue", "Blue"]:
  plts2:= seq
          ( inequal
            ( [ y>psi[j]/eval( f(eta), R)(eta),
                y<psi[j+1]/eval( f(eta), R)(eta)
              ],
              eta=0..5,
              y=-2..2,
              color=cols[j],
              nolines
            ),
            j=1..8
          ):
  plts3:= seq
          ( plot
            ( psi[j]/eval( f(eta), R)(eta),
              eta=0..5,
              color=cols[j]
            ),
            j=1..9
          ):
  display([plts2, plts3]);

Warning, complex or non-numeric value encountered; trying to find a feasible point

Warning, limiting number of function evaluations (1000000) is reached; set initial point equal to extremum point obtained, increase evaluationlimit option and continue search

#
# Generate constraints forcing all variables to be
# in the range 0..1
#
  constr:= [ seq
            ( j=0..1,
              j in `union`
                   ( indets~
                     ( [ seq
                         ( eval(cat(e,j)),
                           j=1..26
                         )
                       ],
                      'name'
                     )[]
                   )
            )
          ];

#
# Solve the problem again, this time with constraints on
# all variables. "Solution" is found much quicker, but
#
# 1. residuals are much higher (not good)
# 2. many of the variables are on the constraint limits -
#    never a good sign
#
  ans := DirectSearch:-SolveEquations
         ( [seq( cat(e,j), j=1..26)],
           constr,
           evaluationlimit = 1000000
         );