@trace I would begin by making m concrete. As many examples as possible.

@trace I would begin by making m concrete. As many examples as possible.

@trace Near the bottom of p.18 the authors make a change of variable. So that has to be done.
Here are the beginnings for the case considered m(r) = n-1.
I start as before.
restart;
f := (x, y, z)-> -1-z-3*y+x^2-x*z;
g := (x, y, z)-> x*z/m(z/y)-2*z*y+4*y+x*y;
h := (x, y, z)-> -x*z/m(z/y)-2*z^2+4*z;
m:=s->n-1;
res := solve({f(x, y, z) = 0, g(x, y, z) = 0, h(x, y, z) = 0}, {x, y, z});
pts:=map(subs,[res],[x,y,z]);
J := unapply(VectorCalculus:-Jacobian([f(x, y, z), g(x, y, z), h(x, y, z)], [x, y, z]), x, y, z):
use ev=LinearAlgebra:-Eigenvalues in
  for p in pts do
  ev(J(op(p)))
  end do
end use;
#For convenience I'm using x1,x2,x3 instead of x,y,z:
sys := diff(x[1](t), t) = f(x[1](t),x[2](t), x[3](t)), diff(x[2](t), t) = g(x[1](t),x[2](t), x[3](t)), diff(x[3](t), t) = h(x[1](t),x[2](t), x[3](t));
p:='p':
alias(P=add(p[i](t)^2,i=1..3)); #To make output shorter
varchange:={seq(x[j](t)=p[j](t)/(1-sqrt(add(p[i](t)^2,i=1..3))),j=1..3)};
PDEtools:-dchange(varchange,{sys},[seq(p[i](t),i=1..3)]);
solve(%,{seq(diff(p[i](t),t),i=1..3)});
psys:=collect(%,[p[1](t),p[2](t),p[3](t)],distributed,factor);
#This is the system in Poincare coordinates. You need to specify n, solve numerically and show plots in the (p1,p2)-plane.

@Markiyan Hirnyk The answer returned by fsolve doesn't mean anything here. Try

evalf[40](eval(J,{z=0.1e-2,x=-12.30640295}));
evalf[40](eval(J,{z=0.1e-2,x=-12.4}));

The problem surely is that J itself is extremely small at those x-values. The "roots" are not roots at all.

My argument was and remains: Since diff(J,x) < 0 for x<0 and limit(J,x=-infinity) = 0 the equation J = 0 has no negative root for x for any given value of z. 

@Markiyan Hirnyk The answer returned by fsolve doesn't mean anything here. Try

evalf[40](eval(J,{z=0.1e-2,x=-12.30640295}));
evalf[40](eval(J,{z=0.1e-2,x=-12.4}));

The problem surely is that J itself is extremely small at those x-values. The "roots" are not roots at all.

My argument was and remains: Since diff(J,x) < 0 for x<0 and limit(J,x=-infinity) = 0 the equation J = 0 has no negative root for x for any given value of z. 

@orawajar You can use the Routh-Hurwitz criterion. For a polynomial with real coefficients this gives necessary and sufficient conditions for all roots to have negative real parts.
In case of a polynomial of 4th degree the conditions are:

p:=LinearAlgebra:-CharacteristicPolynomial(A,lambda);
for i from 0 to 4 do a[i]:=coeff(p,lambda,i) end do;
All the coefficients have to be positive and furthermore
a[3]*a[2]>a[4]*a[1];
a[3]*a[2]*a[1]>a[4]*a[1]^2+a[3]^2*a[0];
#Not surprising then that the following returns FAIL:
PolynomialTools:-Hurwitz(p,lambda);

@orawajar You can use the Routh-Hurwitz criterion. For a polynomial with real coefficients this gives necessary and sufficient conditions for all roots to have negative real parts.
In case of a polynomial of 4th degree the conditions are:

p:=LinearAlgebra:-CharacteristicPolynomial(A,lambda);
for i from 0 to 4 do a[i]:=coeff(p,lambda,i) end do;
All the coefficients have to be positive and furthermore
a[3]*a[2]>a[4]*a[1];
a[3]*a[2]*a[1]>a[4]*a[1]^2+a[3]^2*a[0];
#Not surprising then that the following returns FAIL:
PolynomialTools:-Hurwitz(p,lambda);

@Axel Vogt Right, if the normal/expanded line wasn't there things work out fine in this case.
However, with that line:

eval(1/(2^s-1), s = x - 2*Pi*I*3/log(2));
normal(%,expanded);
series(%, x=0);
coeff(%, x, -1);
Result 0.

@Alejandro Jakubi I forgot to take care of the case s = infinity. Translating to zero makes the code somewhat simpler, but not much. Set x,x0:=op(a) and after assigning to g in the infinity case set x0:=0, then make the trivial adjustments.

@laetititia Heaviside(t) = 1 for t > 0 and  Heaviside(t) = 0 for t < 0. In Maple it is left undefined for t = 0, but the user can change that.
Every piecewise defined function can be expressed in terms of translated versions of Heaviside.
In your example:
p:=x->piecewise(x<10,2,x<11,10,2);
convert(p(x),Heaviside);
#Result: -8*Heaviside(-11+x)+8*Heaviside(-10+x)+2
#You can convert this back tp piecewise:
convert(%,piecewise);
#Here you will see some undefined values due to H(0) not being defined.

It looks like a command line plot, but how you got it in the Standard interface I don't know. I have seen this in the past very rarely, but never found out what caused it.
Have you tried restart? Or if that doesn't make the problem go away, try closing Maple and opening it again.

@Alejandro Jakubi Yes, a similar problem also raised by Marko Riedel is found here:

http://www.mapleprimes.com/questions/141286-Computing-A-Residue

Another choice is to change residue only:

Residue:=subs(series=MultiSeries:-series,eval(residue));

or (in this case)
Residue:=subs(normal=simplify@normal,eval(residue));

One more:

Why make the variable change in residue from say s = s0 to x = s-s0 as is done in the code?
Why not just like this:
residue:=proc (f, a::(name = anything)) local g, i, t, x, t1;
   options remember, `Copyright (c) 1990 by the University of Waterloo. All rights reserved.`;
   x := op(1, a);
   if has(op(2, a), x) then
      error "invalid point"
   elif type(op(2, a), 'infinity') then
      g := -(eval(f, x = 1/x))/x^2
   else
      g := f #Changed
   end if;
   if hastype(g, 'RootOf') then g := convert(g, 'radical') end if;
   for i from 0 to 5 do
      try
         t := series(g, a, i^2+2) #Changed
         catch:
         t := FAIL
      end try;
      if t = 0 then return 0 end if;
      if t = FAIL or not type(t, 'laurent') then next end if;
      t1 := order(t);
      if t1 = infinity or -1 < t1 then return coeff(t, x-op(2,a), -1) end if #Change in coeff
   end do;
   'residue(args)'
end proc;

@Alejandro Jakubi I think it is a problem with series, not coeff.
Debugging residue we see that g below appears. g is correct, but the series t for g is wrong:

restart;
g := 1/(2^s*2^((6*I)*Pi/ln(2))-1); #From de
g1:=simplify(g);
t:=series(g,s=0,2); #Wrong
series(g1,s=0,2); #Correct
coeff(%,s,-1);
#MultiSeries:-series does a better job:
MultiSeries:-series(g,s=0,2);
coeff(%,s,-1);
simplify(%);

Converting piecewise to Heaviside helps on the first order equation in u(x) = z'(x).
Since a graph is desired we need initial conditions. Here I take z(0) = 0, z'(0) = 0.

restart;
Eq2 := diff(z(x), x, x) = p(x)*sqrt(1+diff(z(x), x)^2);
p:=x->piecewise(x<10,2,x<11,10,2);
#The following doesn't work with or without conversion
dsolve({convert(Eq2,Heaviside),z(0)=0,D(z)(0)=0}); #error
#The first order equation in u(x) = z'(x):
Eq1:=subs(diff(z(x),x)=u(x),Eq2);
dsolve({convert(Eq1,Heaviside),u(0)=0}); #z'(0) = 0
res:=int(rhs(%),x=0..x)+C;
eval(res,x=0);
sol:=eval(res,C=0); #Since z(0) = 0 is required
plot(ln(sol),x=0..20);
#This compares well with the corresponding numerical solution:
res:=dsolve({Eq2,z(0)=0,D(z)(0)=0},numeric);
plots:-odeplot(res,[x,ln(z(x))],0..20);

Note 1: Using simplify/piecewise as Carl does is nicer I think than convert/Heaviside.
Note 2: However, for some reason the Heaviside version handles my initial conditions z(0)=0,z'(0)=0 better than does the simplify/piecewise version:
Eq1a:=simplify(Eq1,piecewise);
dsolve({Eq1a,u(0)=0}); #Two solutions: One of them wrong
combine(%);
zz:=op(indets(%,`local`));
about(zz);
#This could be handled artificially like this:
eval(dsolve({Eq1a,u(0)=epsilon}),epsilon=0);

@Carl Love Thanks for the warning. Even the title (question) is long!