You have

pde33 := subs({D(f)(w)*(D^(3)(f)(w) = -c/3*(D^(4)(f)(w)}, pde3);

Those D^3 and D^4 should be D@@3 and D@@4.

After fixing a couple of misplaced parenthesis, the correct version of the above becomes

pde33 := subs({D(f)(w)*(D@@3)(f)(w) = -c/3*(D@@4)(f)(w)}, pde3);

The first version of the code that you posted is pretty much complete but it is carelessly written.  In Maple, there is a major difference between the "=" and ":=" operators.  The purpose of the first one is to state the equality of the left-hand and right-hand operands.  The second one is to define a symbol in terms of other symbols.   That's the first concept that you need to learn and master before you begin making use of Maple.

In the code that you have presented, those two operators are used randomly throughout, rendering the code meaningless.  For instance, where you have

x = a + i*h;

what you really mean is to define x to be a short-hand notation for a + i*h.  Therefore, that should be

x := a + i*h;

Your first task then is to go through that code very carefully and figure out where you need "=" and where ":=".

After you have fixed that, you need to change the few occurrences of rem(j,2).  What you need there is irem(j,2).  

If you do those two steps correctly, you will have a working proc.  Next, you need to learn how to use a proc.

Your proc is designed to integrate a function f(x,y) on a domain the x-y plane bounded by the lines x=a and x=b, and the curves y = c(x) and y=d(x).  To use the proc, you need to supply it with the three functions f, c, and d.  For instance, if f(x,y) = cos(y), c(x) = 0, and d(x) = x, you need

f := (x,y) -> cos(y);
c := x -> 0;
d := x -> x;
APV(f, 0.0, 1.0*Pi, c, d, 6, 6);

If you have been careful, this will produce the answer 2.000980454 which is pretty close to the exact solution which is 2.

I am not posting my corrected version of the code because I don't want to deprive you of the benefits of learning from this exercise.

You have α=0.1 and r goes from −1 to 1. What, in your opinion, is the value of r^(alpha−2) when r=−0.5?

I have fixed quite a few errors in your worksheet, most important of which are:

1. The ":=" and "=" mean quite different things to Maple.  The mathematical equal sign is "=".
2. The derivative at time t=0 should be written as D(xi)(0).

The numbers in your equations are quite odd.  Do you really want to deal with functions such as cos(10^12*t)?  That function oscillates billions of times in the time interval 0 < t < 1.

Anyway, here is the corrected worksheet, for whatever it's worth.
 

Download 66-corrected.mw

restart;
with(plots):
de := dy/dx = -y^2*(5*y-1)^2;
solve(rhs(de));
solve(rhs(de)<0);
arrowhead := proc(head::list)
	local x, y, L;
	x, y, L := op(head);
	plottools:-polygon([[x, y-L/3], [x+L,y], [x,y+L/3], [x+L/6,y]], _rest);
end proc:
arrowheads := proc(heads::listlist)
	seq(arrowhead(head, _rest), head in heads);
	display(%);
end proc:
display(
	pointplot([[-0.25,0],[0.4,0]], connect),
	pointplot([[0,0], [1/5,0]], symbol=solidcircle, symbolsize=35),
	arrowheads([[-0.15,0,-0.035], [ 0.1,0,-0.035], [ 0.35,0,-0.035]], color=red),
	textplot([[0,-0.05, typeset(0)], [1/5, -0.05, typeset(1/5)]], font=[times,bold,16]),
NULL, scaling=constrained, axes=none, size=[500,150]);

Download worksheet: mw.mw

The error message is due to the misplaced argument 'maxerror' in your call to minimax.  The fourth argument of minimax is expected to be a prescribed weight function.  The fifth argument can be 'maxerror' if you want.

Letting the weight function to be 1, we get

minimax(Z(x), x = 1.666666666667 .. 3.5, [2, 1], 1, 'maxerror')

This leads to the error message

Error, (in numapprox:-remez) error curve fails to oscillate sufficiently; try different degrees

So let's change the degree:

minimax(Z(x), x = 1.666666666667 .. 3.5, [3,1], 1, 'maxerror')

This produces the desired result.

Or try with the default weight, which is Z(x):

minimax(Z(x), x = 1.666666666667 .. 3.5, [3,1], Z(x), 'maxerror')

whch also works and produces a smaller maxerror.

PS: Your original message seems to indicate that your ultimate goal is to calculate a certain integral, and that you are attempting to approximate your piecewise function with a rational function in order to evaluate that integral.  That doesn't look like a good approach to me.  If you state the integral that you wish to calculate, someone here may suggest a more straightforward approach.

Why not solve the two equations together as a system?

[I promoted this from a reply to an answer since it is the optimum solution IMO - @dharr ]

Your PDE has two independent variables, x and t. In your formulation, you are changing x but not t, but PDEtools:-dchange expects that you change both independent variables, let's say from (x,t) to (xi,tau).  You want tau = t which is okay; you just have to tell dchange that that's what you want. Therefore we change your code to:

restart;
tr := { x = xi + mu*tau, t = tau, u(x,t) = U(xi) };
pde := diff(u(x,t),t) +p*u(x,t)*diff(u(x,t),x) + q* diff(u(x,t),x$3) = 0;
PDEtools:-dchange(tr, pde, [xi, tau, U]);

This yields the expected result

restart;
with(plots):
arrowhead := proc(head::list)
	local x, y, L;
	x, y, L := op(head);
	plottools:-polygon([[x, y-L/3], [x+L,y], [x,y+L/3], [x+L/6,y]], _rest);
end proc:
arrowheads := proc(heads::listlist)
	seq(arrowhead(head, _rest), head in heads);
	display(%);
end proc:
display(
	plot(x^2/4, x=-4..4, thickness=3, color=blue),
	plot([0,lambda, lambda=-2..4], thickness=3, color=blue),
	plot(
		[[[-4,-1],[4,-1]], [[-4,0],[4,0]], [[-4,1],[4,1]], [[-4,2],[4,2]]], 
		linestyle=dash, color="Green"),
	pointplot([[0,-1], [0,0], [-2,1], [0,1], [2,1],
		[-2*sqrt(2),2], [0,2], [2*sqrt(2),2]],
		symbol=solidcircle, symbolsize=20, color=black),
	arrowheads([
		[-2,-1,-0.3], [2,-1,-0.3],
		[-2.7,0,-0.3], [2.7,0,-0.3],
		[-3.3,1,-0.3], [3.3,1,-0.3],
		[-1.0,1,0.3], [1.0,1,0.3],
		[-3.8,2,-0.3], [3.8,2,-0.3],
		[-1.5,2,0.3], [1.5,2,0.3]
	], color=red),
axes=boxed, labels=[x,lambda], view=-2..3, scaling=constrained);

Download worksheet: mw.mw

restart;
with(plots):
eq1 := 2*x^2 + 14*y^2 = 7;
eq2 := 2*x^2 - 2*y^2 = 3;

#Intersection point in the first quadrant
pt := fsolve({eq1,eq2}, {x=0..2, y=0..1});

# Calculate the normal vector to eq1, then turn it by 90 degrees
# to make it into a tangent vector.
# Additionally, normalize the tangent vector to unit length.
n1 := < diff(lhs(eq1),x), diff(lhs(eq1),y) >:
< -%[2], %[1] >:
t1 := % / sqrt(%^+ . %);

# Repeat the calculations for eq2:
n2 := < diff(lhs(eq2),x), diff(lhs(eq2),y) >:
< -%[2], %[1] >:
t2 := % / sqrt(%^+ . %);

display(
	implicitplot([eq1, eq2], x=-2..2, y=-1..1, color=[red,blue]),
	arrow(eval([<x,y>, 0.8*t1], pt)[], color=red),
	arrow(eval([<x,y>, 0.8*t2], pt)[], color=blue),
	pointplot(eval(<x,y>, pt), symbol=solidcircle, symbolsize=24, color=white),
	pointplot(eval(<x,y>, pt), symbol=circle, symbolsize=24, color=black),
	scaling=constrained);

Download worksheet: mw.mw

Change
S1 := plot3d(F, x = -Pi/2 .. Pi/2, t = 0 .. 2*Pi, color = "Green");
to
S1 := plot3d(F, x = -Pi/2 .. Pi/2, t = 0 .. 2*Pi, color = "Green", transparency=0.4);

Change
S2 := plot3d(G, x = 0 .. 100, t = 0 .. 2*Pi, color = "Cyan");
to
S2 := plot3d(G, x = 0 .. 1, t = 0 .. 2*Pi, color = "Cyan");

I assume that the first r2 in your equation is meant to be r1.

If so, then the data is inconsistent since for all x and y, the left-hand side of the equation is always greater than the right-hand side, and therefore the equality cannot hold.

You can see that in the picture below where I plot the equation's left- and right-hand sides for a range and x and y.

restart;
eq := (m1 + m2)*(x^2+y^2) + 2*m1/r1 + 2*m2/r2 = c;
EQ := subs(
  r1 = sqrt((x-d1)^2 + y^2 + z^2),
  r2 = sqrt((x-d2)^2 + y^2 + z^2),
  d1 = 0.6, d2 = 0.4,
  m1 = 10, m2 = 2, c = 20,
  z = 0,
eq);
plot3d([rhs(EQ),lhs(EQ)], x=-1..2, y=-1..1,
    numpoints=1000, view=0..300, color=[gold, cyan], style=patchcontour);

Depending on the way you want to use the indices, this alternative to Carl Love's suggestion may also be useful:

op(node[2,3,5]);

which yields 2, 3, 5.

It is possible to plot those curves through solving differential equations.but that's too much work for little gain.  Maple's contourplot function is much more effective for that purpose.  Here is how to produce the curves at the bottom of that website:

restart;
with(plots):
# Some judicious selection of contour levels:
levels1 := [ i^2 $i=0..6 ] /~ 5 +~ 1;
levels2 := i/10 $i=1..9;
levels := [levels1[], levels2[]];
display(
    contourplot(exp(y^2)*cos(x)^2, x=-Pi..Pi, y=-2..2, contours=levels),
    contourplot(y*sin(x), x=-Pi..Pi, y=-2..2, contours=11,
        filledregions=true, coloring=[red,yellow]),
    scaling=constrained, axes=none);

Download worksheet: mw.mw

restart;
with(plots):
with(LinearAlgebra):
f := x -> x^2/2 - 1;
g := x -> 2 - x^2/4;

F := <f(x), x*cos(t), x*sin(t) >;

S1 := plot3d(F, x=-2..2, t=0..2*Pi, color="Orange"):

F, 1.5*Normalize(diff(F,x), 2):
eval(%, {x=2,t=Pi}):
A1 := arrow(%, color=red):

F, 1.5*Normalize(diff(F,t), 2):
eval(%, {x=2,t=Pi}):
A2 := arrow(%, color=red):

G := < g(x), x*cos(t), x*sin(t) >;

S2 := plot3d(G, x=-2..2, t=0..2*Pi, color="Cyan"):

G, 1.5*Normalize(diff(G,x), 2):
eval(%, {x=2,t=Pi}):
A3 := arrow(%, color=red):

eval(F, t=0):
convert(%, list):
T1 := tubeplot(%, x=-2..2, radius=0.05, style=surface, color=yellow):

eval(G, t=0):
convert(%, list):
T2 := tubeplot(%, x=-2..2, radius=0.05, style=surface, color=blue):

display(S1, S2, A1, A2, A3, T1, T2, scaling=constrained,
        axes=framed, labels=[x,y,z]);