There are two ways to handle this. The first way involves checking mathematically whether the modular inverse exists before trying to compute it. The second way is to trap the error.

The mathematical way: For integers a, b with b > 1, the modular inverse 1/a mod b exists iff a and b are relatively prime. In Maple, we check that igcd(a,b) = 1. I'll assume that you never have a ridiculous situation such as b=0 or b=1; so I won't check for those.

for i to s1 do
     for c2 to 15 do
          if
               igcd(H+1, A[i]) = 1 and
               `mod`(c2+(J+H)*(`mod`(1/(H+1), A[i])), A[i]) = 0
          then
               d[c2]:= d[c2]+l(A[i])
          end if
     end do
end do;

If the first part of an and clause evaluates to false, then the second part is never evaluated or executed. That's why the above code is safe.

The lazy error-trapping way: This requires using the Maple statement try/catch.

for i to s1 do
     for c2 to 15 do
          try
               if `mod`(c2+(J+H)*(`mod`(1/(H+1), A[i])), A[i]) = 0 then
                    d[c2]:= d[c2]+l(A[i])
               end if
          catch "the modular inverse":
               #Ignore the error; don't change d[c2].
          end try
     end do
end do;

In both cases I've assumed that you intended l(A[i]) to represent a function (or procedure) call to a function named l. If you actually intended multiplication then you need to make that l*A[i].

Note that p/q mod b may exist even though 1/q mod a does not. In light of this you may wish to extend the domain of your computation (and simplify the code) thus:

for i to s1 do
     for c2 to 15 do
          try
               if c2+(J+H)/(H+1) mod A[i] = 0 then
                    d[c2]:= d[c2]+l(A[i])
               end if
          catch "the modular inverse":
               #Ignore the error; don't change d[c2].
          end try
     end do
end do;

This is an oft-reported bug with curves in 3d plots that only affects some users. A workaround is to include linestyle= solid in the plot command.

plots:-contourplot3d(-5*x/(x^2+y^2+1), x= -3..3, y= -3..3, filledregions, linestyle= solid);

plot3d(-5*x/(x^2+y^2+1), x= -3..3, y= -3..3, style= contour, contours= 6, linestyle= solid);

assign(seq(seq(psi[4*(i-1)+j+1]= psi[i,j], i= 1..2), j= 0..3));

type(z1, specfunc(anything, `.`));

                               true

or

type(z1, ('`.`')(anything, anything));

                               true

You should clearly state the question. Is it to maximize the profit (or minimize the loss)? Is it to find the break-even point?

I think that there's missing information from your problem, or that you are mixing units incorrectly (such as dollars and cents). I think that because, as stated, the problem does not have a positive real solution.

price:= x-> 50-x:
Cost:= x-> 120+30*x:
Profit:= x-> x*price(x) - Cost(x):
maximize(Profit(x), location);

So the minimal loss is $20 at a production of 10 units and a price of 50 - 10 = $40.

It is possible, with algcurves[parameterization], to get a single parameterization instead of four.

restart:
f1:= 3-2*x^2-y^2:  f2:= x^2+2*y^2:
(X,Y):= algcurves[parametrization](f1-f2, x, y, t)[];

Z:= eval(f1, [x,y]=~ [X,Y]):
evalf(Int(sqrt(diff(X,t)^2+diff(Y,t)^2+diff(Z,t)^2), t= -infinity..infinity));

                7.640395578

We can get the exact value of the integral via a trigonometric substitution. The integrand simplfies remarkably.

S:= simplify(eval([X,Y,Z], t= tan(s)));

Int(sqrt(simplify(add(diff(S[k],s)^2, k= 1..3))), s= -Pi/2..Pi/2);
value(%);

evalf(%);

                7.640395576

plots:-spacecurve(S, s= -Pi/2..Pi/2, linestyle= solid, axes= frame, thickness= 3);

Note that the coefficients of both derivatives are 0 at x=0, which is thus a singular point. Thus it is not proper to give initial conditions at that point.

Here is a way to automate this computation somewhat with Maple. For this to work it is important to use a consistent ordering of the letters in the variable names.

n, nC, nH, nB, nHB, nCB, nCH, nCHB:= 800, 224, 240, 336, 64, 80, 40, 24:
eval(add(cat('n', s[])*(-1)^irem(nops(s),2), s= combinat:-powerset([C,H,B])));
                      160

A:= <10,15>:
B:= <12,8>:
Req:= <100,100>:
LinearAlgebra:-LinearSolve(< A | B >, Req);

Here is a procedure to count the minimal number of transpositions to achieve a permutation. It does not error-check the input, i.e., it does not check that the first list has unique elements or that the second list is a permutation of the first.

Transpositions:= proc(L1::list, L2::list)
description "Count the transpositions in a permutation.";
local x, y, pos:= table(L1 =~ [$1..nops(L1)]);
     add(nops(x)-1, x= convert([seq(pos[y], y= L2)], disjcyc))
end proc:

Example of use:

Transpositions([w,x,y,z], [y,x,z,w]);

                   2

Here's my Maple implementation of the pseudo-code that you gave:

GreedyVertexColoring:= proc(G::GraphTheory:-Graph)
uses GT= GraphTheory;
local
     #Initialize colors to 0.
     Colors:= proc(v) option remember; 0 end proc,
     N, v, k
;
     for v in GT:-Vertices(G) do
          N:= map(Colors, {GT:-Neighbors(G,v)[]});
          for k while _C(k) in N do end do;
          Colors(v):= _C(k)
     end do;
     # Return remember table.
     op(4, eval(Colors))
end proc:

Example of use:

G:= GraphTheory:-SpecialGraphs:-PetersenGraph():
GreedyVertexColoring(G);

Colors:= table([_C(1)=red, _C(2)=green, _C(3)=blue]):
for v in GraphTheory:-Vertices(G) do
     GraphTheory:-HighlightVertex(G, v, Colors[C[v]])
end do:
GraphTheory:-DrawGraph(G);

I'll do this one because there are some idiosyncratic-to-Maple stumbling blocks to getting the proof.

Download sequence.mw

You almost have it. Your problem is that there is no need for any looping command. Some local variables will help also.

First_Principle:= proc(f::appliable)
local x, h, A:= simplify((f(x+h)-f(x))/h);
     Limit(A, h = 0) = limit(A, h = 0)
end proc:

Example of use:

First_Principle(sin);

You can remove from the solutions those for which the evaluation of R leads to error like this:

Error:= proc() error NumericException(division_by_zero) end proc:
NumericEventHandler(division_by_zero= Error):
CheckSolution:= proc(sol::set(name=anything))
     try
          eval(R, sol)
     catch "numeric exception: division":
          return true
     end try;
     false
end proc:

remove(CheckSolution, sols);

Nothing special is required because of the change of version. Just open it as you would open any other file.