@mapleatha  In Maple 2015-2017 everything works properly. Execute  restart  command first:

restart;
subs([A=1/6,B=8/15,C=3/10], A/(s+1)+B/(s-2)+C/(s+3));

@max125 

1. Try yourself to find errors in your procedures.

2. Of course, this procedure can be written much shorter, for example:

Orthocenter:=proc(A::list, B::list, C::list)
local H;
if (B[1]-A[1])*(C[2]-A[2])=(B[2]-A[2])*(C[1]-A[1]) then error "The points A, B, C lie on one straight line" fi;
H:=proc(A, B, C)
local N;
N:=B-A;
N[1]*('x'-C[1])+N[2]*('y'-C[2])=0;
end proc;
solve({H(A,B,C), H(B,C,A)});
end proc:

Example of use:

Orthocenter([-3,3], [1,3], [-3,0]);

                                                  {x = -3, y = 3}

The problem is to symbolically prove that the determinant is 0. This is trivial:

LinearAlgebra:-Determinant(M):
is(%=0);
                                                              true

Another option:

simplify(LinearAlgebra:-Determinant(M));
                                                                0

@digerdiga  indexed  is the type of, for example,  a[n] :

whattype(a[n]);
                                  indexed

op(a[n]);
                                       n

The  evalindets  command  (see above)  applies  the procedure  t->op(t)*p+q+b[op(t)]  to any operand  t  of  eq1  of the  indexed  type. 

@Thomas Richard  Thank you very much! The third item on the list "Questions" helped.

I have the same problem in 32 bit Maple 2017.3 on Windows 10. The problem appeared about 2.5 months ago. Since then, I have reinstalled Maple and even been able to update it to Maple 2017.3, but the problem has remained (the standard GUI hangs on boot).

@pgl771  This is even simplier:

plot3d(y^2+x, x = 0 .. 2, y = x-2 .. x, orientation = [220, 75]);

@Les  For plotting in geom3d  package the  draw  command should be used:

draw(pp);

@Carl Love  The Iterator package was introduced in Maple 2016, therefore, unfortunately, OP will not be able to check the operation of your procedure.

@mapleatha  I was able to get only this option (with lambda[1]  in front of I):

restart;
with(InertForm):
Display(`%+`(A,`.`(lambda[1], convert(I,symbol))))=<a,b; c,d>;

@nunavutpanther   You can write this shorter:

L1:=[y,x*y];
L2:=[x^2,x*y];
{L1[], L2[]};

@alfarunner  Do not forget  restart  before running the code. For me in Maple 2015 - 2017 everything works properly:

restart;
z:=(-1/(4*(-u+L)))*(4*p^2*b^2*d-3*g*c*p^2);
u:=2:  L:=1:  c:=0.6: 
constr:={ 1-(2*d)/(p*(sqrt(g)+1))<=b, b<=1, 0<=d, d<=L, c<=p, p<=u, 1<=g, g<=4};  
Optimization:-Maximize(z, constr);

I tried to apply your procedures (in Maple 2017.3 and 2016.2) to the same example and got a strange error:

restart;
ranges:=proc(simpledom::list(relation), X::list(name))
local rez:=NULL, r,z,k,r1,r2;
if nops(simpledom)<>2*nops(X) then error "Domain not simple!" fi;
for k to nops(X) do    r1,r2:=simpledom[2*k-1..2*k][]; z:=X[k];
  if   rhs(r1)=z and lhs(r2)=z then rez:=z=lhs(r1)..rhs(r2),rez; #a<z,z<b
  elif lhs(r1)=z and rhs(r2)=z then rez:=z=lhs(r2)..rhs(r1),rez  #z<b,a<z
  else error "Strange order in a simple domain" fi
od;
rez
end proc:

MultiIntPoly:=proc(f, rels::list(relation(ratpoly)), X::list(name))
local r,rez,sol,irr,wirr, rels1, w;
irr:=[indets(rels,{function,realcons^realcons})[]];
wirr:=[seq(w[i],i=1..nops(irr))];
rels1:=eval(rels, irr=~wirr);
sol:=SolveTools:-SemiAlgebraic(rels1,X,parameters=wirr):
sol:=remove(hastype, eval(sol,wirr=~irr), `=`); 
add(Int(f,ranges(r,X)),r=sol)
end proc:

R := [ (x-4)^2+y^2 <= 25, x^2+(y-3)^2 >= 9, (x+sqrt(7))^2+y^2 >= 16 ];
MultiIntPoly(1, R, [x,y]): 
radnormal(simplify(value(%)));

@mapleatha  In this case use  Standard GUI  in which there are no any parentheses around the exponent.

In Classic interface, you can do so:

expa := 2: 
exp(cat(expa,` t`));

Edit.

@vv  For better perception, I did everything on one plot, eliminated some repetitions in the code and added the  fieldstrength = log   option that equalizes the lengths of the arrows:

restart;
with(plots): with(plottools):
poly:=x>=0, y<=1, y>=-1, y<=4-4*x, y>=4*x-4:
F,G:=piecewise~(`and`(poly), [x^2+y^2,3*sin(x)-y^2], undefined)[]:
display(polygon([[0, -1], [3/4, -1], [1, 0], [3/4, 1], [0, 1]], style = line, thickness = 2), fieldplot([F,G], x =0 .. 1, y = -1 .. 1, arrows = SLIM, color = x, fieldstrength = log));