restart

a := 1;

h := 1;

# pA := [0, 0, 0]; # useless

pB := [a, 0, 0];

pC := [a, a, 0];

pD := [0, a, 0];

pS := [0, 0, h];

local D:
symbols := [B, C, D, S]:

use geom3d in
  seq(point(s, p||s), s in symbols);
  plane(SBC, [S, B, C]);
  plane(SCD, [S, C, D]);
  sbc    := < NormalVector(SBC) >;
  scd    := < NormalVector(SCD) >;

  # from the reference you give :
  # "the dihedral angle, φ... satisfies 0 ≤ φ ≤ π/2"
  # So 2*π/3 and 3*π/4 cannot be solutions as you say.
  
  cosphi := abs(Statistics:-Correlation(sbc, scd));
end use:

phi := identify(arccos(cosphi))

restart

eq_7 := diff(G(eta), eta$2) + sigma*G(eta) = nu

negsol := rhs( ( dsolve(eq_7) assuming sigma < 0) ):
negsol := expand(convert(negsol, trig)):


CS      := [C, S]:
replace := convert(indets(negsol, function), list) =~ CS:
negsol  := collect(eval(negsol, replace), CS):
B1B2    := solve({coeff(negsol, C)=B1, coeff(negsol, S)=B2}, [_C1, _C2]):

eq_8 := eval(eval(negsol, %[]), (rhs=lhs)~(replace));

use IntegrationTools in

  # multiply both member of eq_7 by G'

  expand~(eq_7*~diff(G(eta), eta));
 
  # Integrate over eta

  map((Expand@Int), %, eta);
 
  # evaluate

  value~(%);
 
  # divide each term by G^2 (assuming G <> 0)

  expand(% /~ G(eta)^2);
 
  # isolate (G'/G)^2

  ref := isolate(%, (diff(G(eta), eta) / G(eta))^2);
 
  # add an integration "constant"

  ref := lhs(ref) = rhs(ref) + __K;
 
  # plug in ref the expression of G(eta) given by eq_8

  aux := eval(ref, G = (eta -> eq_8));

  # solve for __K

  K := simplify(solve(aux, __K));

  # Identify the denominator of K

  den := simplify(denom(K), {eq_8=G(eta)});

  # Rewrite K according to theprevious result

  K := numer(K)/den;

  # Finally

  eq_9 := eval(ref, __K=K);

end use:

eq_9;

eq_10 := eval( rhs( ( dsolve(eq_7) assuming sigma > 0) ), [_C1=B1, _C2=B2])

use IntegrationTools in
 
  # plug in ref the expression of G(eta) given by eq_8

  aux := eval(ref, G = (eta -> eq_10));

  # solve for __K

  K := simplify(solve(aux, __K), size):
  K := simplify(numer(K))/denom(K);

  # Identify the denominator of K

  den := simplify(denom(K), {eq_10=G(eta)});

  # Rewrite K according to theprevious result

  K := numer(K)/den;

  # Finally

  eq_11 := eval(ref, __K=K);

end use:

eq_11

zerosol := dsolve(eval(eq_7, sigma = 0)):
eq_12   := rhs(eval(zerosol, {_C1=B1, _C2=B2}));

use IntegrationTools in

  # multiply both member of eq_7 by G'

  expand~(eval(eq_7, sigma=0)*~diff(G(eta), eta));
 
  # Integrate over eta

  map((Expand@Int), %, eta);
 
  # evaluate

  value~(%);
 
  # divide each term by G^2 (assuming G <> 0)

  expand(% /~ G(eta)^2);
 
  # isolate (G'/G)^2

  ref := isolate(%, (diff(G(eta), eta) / G(eta))^2);
 
  # add an integration "constant"

  ref := lhs(ref) = rhs(ref) + __K;
 
  # plug in ref the expression of G(eta) given by eq_8

  aux := eval(ref, G = (eta -> eq_12));

  # solve for __K

  K := simplify(solve(aux, __K), size):
  K := simplify(numer(K))/denom(K);

  # Identify the denominator of K

  den := simplify(denom(K), {eq_12=G(eta)});

  # Rewrite K according to theprevious result

  K := simplify(numer(K)/den);

  # Finally

  eq_13 := eval(ref, __K=K);

end use:

eq_13

restart

# Access to the help page

?sum

# from ,sum help page, example 4 (Maple 2015)

sum(1/k!,k=0..infinity);
evalf(%);

# A first example invoking the exponential function

sum(exp(-k), k=0..infinity);

# Use of the Inert form ('Sum'):

Sum(exp(-k), k=0..infinity) = sum(exp(-k), k=0..infinity);

# Be careful: functions 'add' and 'sum' do mean the same thing.
# This is a common error.
#
# To understand the difference between them: 'add' stands for 'addition' while
# 'sum' stands for... 'sum'
# (see these terms in Maple's dictionary: within the "Maple xxx Help" window: in the
# left pane click on the drop-down box labeled "Page Types",  unclick "Select all"
# and select only "Dictionary").
#
# Once done type addition in the top-left field and "Search": the right page should
# the definition of the term addition and refer to the function contain 'add'.
#
# Next type sum instead of addition and read carefuly the second item in the right pane:
#    2.  The limit of the sequence  of partial sums of the first n terms of an infinite series ,
#        as n tends to infinity.
#

restart:

with(geometry):

with(plots):

Bl := color = black;

y0 := x -> -ln(1 - exp(-x)):

y1 := x -> -ln(-1 + exp(-x)):

y2 := x -> -ln(1 + exp(-x)):

p  := plot(y0(x), x = 0.02 .. 4, scaling = constrained, color = blue):

p1 := plot(y2(x), x = -4 .. 4, scaling = constrained, color = green):

p2 := plot(y1(x), x = -4 .. 0, scaling = constrained, color = red):

display({p, p1, p2}, view = [-4 .. 4, -4 .. 5]);

# First way

plots:-display(
  plots:-inequal({y >= y2(x), y <= y1(x)}, x=-4..0, y=-4..4),
  plots:-inequal({y >= y2(x), y <= y0(x)}, x=0..4, y=-4..4)
);

# A less efficient way but having the characteristic function of the blue domain
# might be interesting.

interior := (x, y) -> piecewise(
                        x >=0
                        , piecewise(y < y0(x) and y > y2(x), 1, 0)
                        , piecewise(y < y1(x) and y > y2(x), 1, 0)
                      ):

plots:-contourplot(
  interior(x, y)
  , x=-4..4, y=-4..4
  , contours=1
  , grid=[200, 200]    # use grid=[400, 400] for a better rendering
  , filled=true
)

Area := int(int(1, y=y2(x)..y0(x)), x=0..4) + int(int(1, y=y2(x)..y1(x)), x=-4..0);
evalf(Area);

restart

ddf := proc(F, {n::posint:=1})

  local f, df, idd, idf, rules, i, deg, opi:
  f  := eval(F, indets(F, name) =~ indets(F, name)(__t));
  df := diff(f, __t$n);

  idd, idf := selectremove(has, indets(df), diff);
  idd := convert(idd, list);
  idf := convert(idf minus {__t}, list);

  rules := NULL:
  for i in idd do
    if is(op([0, 0, 0], convert(i, D)) = `@@`) then
      deg := op(2, op([0, 0], convert(i, D)));
      opi := op([0, 1], convert(i, D));
    else
      deg := 1:
      opi := op([1, 0], i);
    end if:
    rules := rules, i = cat(`#mrow(mi("`, opi, `")`, seq(`,mi("&#x27;")`, j=1..deg), `)`);
  end do:
  for i in idf do
    rules := rules, i = cat(`#mi("`, op(0, i), `")`);
  end do;
  
  eval(df, [rules])
end proc:

 

ddf(z^2+y);
ddf(z^2+y, n=1);
ddf(z^2+y, n=2);
ddf(t, n=4);

# Your result: y*z*(2*z*z'+y')*y''+6*y^2*x''''

expr := y * z * ddf(z^2+y, n=1) * ddf(y, n=2) + 6 * ddf(x, n=4) *y^2;