Another way:

x1 := solve(ln(x+2) = -1):   x2 := solve(ln(x+2) = 2):

plots[implicitplot](x*(y-ln(x+2)) >= 0, x = x1 .. x2, y = -1 .. 2, coloring = [blue, white], filledregions, gridrefine = 5, scaling = constrained);

A more automated solution can be get if you use  IntegerPoints  or  IntegerPoints1  procedure from here

eq1 := J+10*A+50*T = 500:

eq2 := J+A+T = 100:

IntegerPoints1({eq1, eq2, A >= 0, J >= 0, T >= 0}, [J, A, T]);

                                               [[60, 39, 1]]

RungeKutta4 procedure solves your problem for equation y'=f(t,y)  with initial condition  y(t0)=y0 .  The procedure was written on the basis of the formulas from wiki . Formal parameters: f  is the function  (t,y)->f(t,y)  defined by the right side of the equation, t0  and  y0  are defined by the initial condition, Tc  is the step,  T is the value of  the variable t  that t[N-1]<T<=t[N] ,  t[N] is the final value of  t , for which the procedure computes  y(t) . The procrdure returnes the list  [[t[0],y[0]], [t[1],y[1]], ... ,  [t[N],y[N]]] .

restart;

RungeKutta4:=proc(f, t0, y0, Tc, T)

local  N, t, y, h, L, n, k1, k2, k3, k4;

N:=ceil((T-t0)/Tc);

t[0]:=evalf(t0);  y[0]:=evalf(y0);  h:=evalf(Tc);  L[1]:=[t[0], y[0]];

for n from 0 to N-1 do

k1:=f(t[n],y[n]);

k2:=f(t[n]+h/2,y[n]+h/2*k1);

k3:=f(t[n]+h/2,y[n]+h/2*k2);

k4:=f(t[n]+h,y[n]+h*k3);

y[n+1]:=y[n]+h/6*(k1+2*k2+2*k3+k4);

t[n+1]:=t[n]+h;

L[n+2]:=[t[n+1],y[n+1]];

od;

convert(L, list);

end proc:

An model example: to solve the equation  y'=sin(t), y(0)=1  (the exact solution  y=2-cos(x) )

RungeKutta4((t, y)->sin(t), 0, 1, 0.1, 5);

plot(%, color=red, thickness=2, scaling=constrained, view=[0..5,0..3], labels=[t,y(t)]);

 Edit. The code has been corrected according to acer's comment.

First, we make the change  g(n)=sqrt(f(n))+1 , and easy to obtain a new recurrence relation

g(n+1)=g(n)+1  with the initial condition g(1)=sqrt(c)+1  (we assume that  f(1)=с ):

restart;

rsolve({g(n+1)=g(n)+1, g(1)=sqrt(c)+1}, g(n));

solve(%=sqrt(f(n))+1, f(n));

factor(subs({sqrt(c)=d,c=d^2},%));

f(n)=subs(d=sqrt(c),%);  # The final result

Carl, unfortunately your approach works only with integer degrees. Here is another solution that works with any degree:

P := 1+sin(e+e^1.5+e^2.5+e^3):

applyrule(e^n::numeric = `if`(n >= 2, 0, e^n), P);

                                        1+sin(e+e^1.5)

plots[implicitplot]((y-4)*(y-x^2) <= 0, x = 0 .. 2, y = 0 .. 4, filledregions = true, scaling = constrained, gridrefine = 5);

Addition 1. For details see help on  plots[implicitplot]  command  with option  filledregions=true 

Addition 2. Another general method you can see  here  (Picture  procedure).

This is only possible if  B  is not yet calculated:

B := Int(A(x), x = -infinity .. infinity);

op(1,B);

Standard way is using of  allvalues  command:

allvalues(SOLL[2]);

I think the initial problem will be easier and more natural to solve by plot3d command with option filled . The  scaling=constrained option is also useful, because it gives an accurate representation of the body proportions. The style = surface  option is not necessary, because in Maple 2015 and 2016 it operates by default:

plot3d([-y^2, x^2], x = 0 .. 1, y = 0 .. x, axes = normal, filled, scaling = constrained, orientation = [-65, 65, 0]);

Syntax

plots:-shadebetween(-y^2, x^2, y = 0 .. x, x = 0 .. 1, axes = frame, filled, style = surface);

plots down a left-handed coordinate system  xyz .

Compare:

A := plots:-shadebetween(-y^2, x^2, y = 0 .. x, x = 0 .. 1, axes = normal, filled):

B := plots:-shadebetween(-y^2, x^2, x = 0 .. 1, y = 0 .. x, axes = normal, filled):

plots[display](<A | B>);

The simple procedure  RandomWalkOnGraph solves the problem:

RandomWalkOnGraph:=proc(G::Graph, s::{symbol,integer}, N::integer)  # s is a vertex, N is the number of vertices passed

local S, i;

uses GraphTheory, RandomTools;

if N=1 then return [s] else

S[1]:=s;

for i from 2 to N do

S[i]:=Generate(choose(Neighbors(G, S[i-1])));

od; fi;

convert(S, list);

end proc:

Example of use:

G := GraphTheory[Graph]( GraphTheory[Trail](1,2,3,4,5,6,4,7,8,2) );

GraphTheory[DrawGraph](G);

RandomWalkOnGraph(G,2,20);

 Edited.

Use polar coordinates:

restart;

f := (r, phi) ->sqrt(25-r^2):

plot3d([r*cos(phi), r*sin(phi), f(r, phi)], r = 3 .. 4, phi = 0 .. 2*Pi, color = khaki, scaling = constrained, view = [-5 .. 5, -5 .. 5, 0 .. 5], axes = normal, labels = [x, y, z]);

Just add the extra line  

local EQU;  

to your procedure:

Z_DAK_FSOLVE := proc( PR , TR ) # This routine solves the Dranchuk - Abou Kassem equation for the  Z  gas compressibility factor using Maple's fsolve routine.

local  EQU;

      RHOR(Z,PR,TR)        := 0.27*PR / (Z*TR) :

      C1(TR)               := 0.3265 - 1.07  /TR - 0.5339/(TR^3) + 0.01569/(TR^4) - 0.05165/(TR^5) :

      C2(TR)               := 0.5475 - 0.7361/TR + 0.1844/(TR^2) :

      C3(TR)               := 0.1056*(-0.7361/TR + 0.1844/(TR^2)):

      C4(RHOR(Z,PR,TR),TR) := 0.6134*(1.+0.721*((RHOR(Z,PR,TR))^2))*(((RHOR(Z,PR,TR))^2)/(TR^3))*exp(-0.721*((RHOR(Z,PR,TR))^2)) :

      EQU                  := Z - ( 1. + (C1(TR))*(RHOR(Z,PR,TR)) + (C2(TR))*((RHOR(Z,PR,TR))^2) - (C3(TR))*((RHOR(Z,PR,TR))^5) + (C4(RHOR(Z,PR,TR),TR)) ) :

      fsolve( EQU , Z) :

      return Z:

end proc :

Example:

A1:=[0, 0, 0]:  A2:=[1, 0, 0]:  A3:=[1, 1, 1]:

P:=plottools[polygon]([A1, A2, A3], color=red):

plots[display](P, axes=normal, labels=[x,y,z], orientation=[10,60]);

 Edited.

Try  ContoursWithLabels  procedure  from  here

l:=nops(block):

for i to l do

block1[i]:=subsop(seq(j=1-block[i][j], j=3..7), block[i]);

od:

block1 := convert(block1, list);

      block1 := [[0, 1, 1, 1, 0, 1, 0, 1], [0, 1, 0, 1, 1, 1, 1, 1], [0, 1, 0, 1, 0, 1, 0, 0], [0, 1, 0, 1, 0, 1, 1, 1], ...]

or even shorter:

block1:=[seq(subsop(seq(j=1-block[i][j], j=3..7), block[i]), i=1..nops(block))];

Edited.