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.

with(LinearAlgebra):

m := RowDimension(Ans):

n := ColumnDimension(Ans):

k := 0:

for i from 2 to m do

if convert([seq(Im(Ans[i, j]) <> 0, j = 1 .. n)], `or`) then k := k+1; L[k] := i  fi;

od:

L := convert(L, list):

DeleteRow(Ans, L);

a := rsolve({a(0) = 2, a(n+1) = a(n)+a(n)^2}, a(n), 'makeproc'):

Example of use:

a(10); 

with(Statistics):

X := RandomVariable(Normal(0, 1));

f := unapply(PDF(X, t), t);

plot(f(t), t = -5 .. 5, color = blue, filled);