#
# Arbitrary values for alpha, c, gamma
#
  pvals:=[alpha=1, c=2, gamma=3];
  poly:=eval( numer(diff(rel, k)), pvals);
#
# Get all roots of poly including complex ones
#
  [fsolve(poly, complex)];
#
# Select only real roots
#
  crit_points:=select( type, %, realcons);

  restart;
#
# Evaluates the argument before generating latex
#
  latex((4*n-1)/9-7/16*n);
#
# Use the inert form of '/'
#
  latex((4*n-1)%/9-7%/16*n);
  

\frac{n}{144}-\frac{1}{9}
\frac{4 n -1}{9}-\frac{7}{16} n

restart;
with(geometry):
with(plots):
_EnvHorizontalName = 'x':
_EnvVerticalName = 'y':
EQ := proc(M, N)
local eq, sol;
eq := simplify(expand((y - M[2])/(x - M[1]) - (N[2] - M[2])/(N[1] - M[1])));
sol := solve(eq, y);
RETURN(y = sol); end proc:
_local(D):
point(A, [-2, 7]):
point(B, [-5, -2]):
point(C, [8, -7]):
point(E, [1, 4]):
EQ([-5, -2], [8, -7]):
point(D, [1, subs(x = 1, rhs(%))]):
dsegment(sgAD, [A, D]):
BD := distance(B, D):
DC := distance(C, D):
triangle(ABC, [A, B, C]):
area(ABC):
triangle(ABD, [A, B, D]):
area(ABD):
triangle(ADC, [A, D, C]):
area(ADC):
is(area(ABD)/area(ADC) = BD/DC):
triangle(EBD, [E, B, D]):
area(EBD):
triangle(EDC, [E, D, C]):
area(EDC):
triangle(AEC, [A, E, C]):
area(AEC):
triangle(ABE, [A, B, E]):
area(ABE):
is(area(ABE)/area(AEC) = BD/DC):
display
( [ draw
    ( [ A(color = black, symbol = solidcircle, symbolsize = 6),
        B(color = black, symbol = solidcircle, symbolsize = 6),
        C(color = black, symbol = solidcircle, symbolsize = 6),
        ABC(color = blue)
      ]
    ),
    textplot
    ( [ [coordinates(A)[], "A"],
        [coordinates(B)[], "B"],
        [coordinates(C)[], "C3"]
      ],
      align = [above, right]
    )
  ],
  axes = none,
  title = "Lemme du Chevron"
);

  restart:
  with(plots):
  Params := [ fw = 0.2, M = 0.5, Q = 0.5, Pr = 6.2, phi = 0.05,
              rf = 997.1, kf = 0.613, cpf = 4179, btf = 0.00003,
              p1 = 0.01, p2 = 0.05, p3 = 0.05, rs1 = 5100, ks1 = 3007.4,
              cps1 = 410, bs1 = 0.0002, rs2 = 2200, ks2 = 5000, cps2 = 790,
              bs2 = 0.0005, rs3 = 3970, ks3 = 40, cps3 = 765, bs3 = 0.0004,
              A1 = B1*p1 + B2*p2 + B3*p3, B1 = 1 + 2.5*phi + 6.2*phi^2,
              B2 = 1 + 13.5*phi + 904.4*phi^2, B3 = 1 + 37.1*phi + 612.6*phi^2,
              B4 = (ks1 + 2*kf - 2*phi*(kf - ks1))/(ks1 + 2*kf + phi*(kf - ks1)),
              B5 = (ks2 + 3.9*kf - 3.9*phi*(kf - ks2))/(ks2 + 3.9*kf + phi*(kf - ks2)),
              B6 = (ks3 + 4.7*kf - 4.7*phi*(kf - ks3))/(ks3 + 4.7*kf + phi*(kf - ks3)),
              A2 = 1 - p1 - p2 - p3 + p1*rs1/rf + p2*rs2/rf + p3*rs3/rf,
              A3 = B4*p1 + B5*p2 + B6*p3, A4 = 1 - p1 - p2 - p3 + p1*rs1*cps1/(rf*cpf) + p2*rs2*cps2/(rf*cpf) + p3*rs3*cps3/(rf*cpf)
          ]:
  ODEs:= [ A1*(diff(f(x), x, x, x))/(A2*phi)-(diff(f(x), x))^2-M^2*(f(x))+f(x)*(diff(f(x), x, x)),
           A4*Pr*phi*(diff(Theta(x), x, x))/A3+f(x)*(diff(Theta(x), x))+Q*Theta(x)
         ]:
  ODEs:=eval[recurse](ODEs, Params):
  BCs := [f(0) = fw, D(f)(0) = 1, Theta(0) = 1, D(f)(1) = 0, Theta(1) = 0]:
  BCs := eval[recurse](BCs,Params):
  sol:=dsolve([ODEs[], BCs[]], numeric);
  odeplot( sol,
           [ [ x, f(x)],
             [x, Theta(x)]
           ],
           x=0..1,
           color=[red, blue]
         );

  restart;
  plot3d( 1,
          theta = 0..Pi,
          phi = 0..2*Pi,
          coords = spherical,
          style =surface,
          scaling=constrained
        );

  restart:
  with(geometry):
  
  Pn := [ (5*x + 12*y)*(-384 + 48*x + 55*y)/(225*x^2 + 225*y^2 - 1800*x - 602*y),
          2*(12*x - 5*y)*(-384 + 48*x + 55*y)/(225*x^2 + 225*y^2 - 1800*x - 602*y)
        ]:
  line(l1, 2+3*x-y, [x,y]):
#
# Define the conic from five points
#
  conic( c,
         [ seq
           ( point
             ( cat(P__, j),
               eval
               ( eval( Pn, isolate( Equation(l1), y ) ),
                 x=j
               )
             ),
             j=1..5
           )
         ],
         [x,y]
        ):
#
# Some properties of the conic
#
  form(c);
  simplify(coordinates(center(c)));
  simplify~(coordinates~(foci(c)));
  simplify(Equation(c));
#
# Draw the conic
#
  c_cen:=coordinates(center(c)):
  vr:= z-> [ c_cen[1]-z..c_cen[1]+z,
             c_cen[2]-z..c_cen[2]+z
           ]:
  draw( c(color=red),
        view=vr(30)
      );

  restart:
  N := 4: M := .1: Kp := .1: Gr := 0.1e-1: Gc := 0.1e-1: Pr := 1: S := 0.1e-1: Sc := .78: Kc := 0.1e-1: La := 1:
  f__N:=  sum((p^(i))*f[i](x), i = 0 .. N) :
  Theta__N:=  sum((p^(i))*Theta[i](x), i = 0 .. N) :
  Phi__N:= sum((p^(i))*Phi [i] (x), i = 0 .. N):
  HPMEq1 := (1-p)*(diff(f__N, x, x, x))+p*(diff(f__N, x, x, x)+(1/2)*(diff(f__N, x, x))*f__N-(M^2+Kp)*(diff(f__N, x)-La)+Gr*Theta__N+Gc*Phi__N):
  HPMEq2 := (1-p)*(diff(Theta__N, x, x))/Pr+p*((diff(Theta__N, x, x))/Pr+(1/2)*(diff(Theta__N, x))*f__N+S*Theta__N):
  HPMEq3 := (1-p)*(diff(Phi__N, x, x))/Sc+p*((diff(Phi__N, x, x))/Sc+(1/2)*(diff(Phi__N, x))*f__N+Kc*Phi__N):

  for i from 0 to N do
      equ[1][i] := coeff(HPMEq1, p, i) = 0:
  end do:
  for i from 0 to N do
      equ[2][i] := coeff(HPMEq2, p, i) = 0:
  end do:
  for i from 0 to N do
      equ[3][i] := coeff(HPMEq3, p, i) = 0:
  end do:
  cond[1][0] := f[0](0) = 0, (D(f[0]))(0) = 0, Theta[0](0) = 1, Phi[0](0) = 1, Theta[0](5) = 0, Phi[0](5) = 0, (D(f[0]))(5) = 1:
  for j to N do
       cond[1][j] := f[j](0) = 0, (D(f[j]))(0) = 0, Theta[j](0) = 0, Phi[j](0) = 0, Theta[j](5) = 0, Phi[j](5) = 0, (D(f[j]))(5) = 0:
  end do:

  ans:=dsolve({equ[1][0], equ[2][0], equ[3][0], cond[1][0]}):
  for k from 1 by 1 to 4 do
      ans:=`union`( ans, dsolve( { eval({equ[1][k], equ[2][k], equ[3][k]}, ans)[], cond[1][k]})):
  od:
  Phi:= evalf(eval(Phi[k-1](x), ans));
  Theta:= evalf(eval(Theta[k-1](x), ans));
  f:= evalf(eval(f[k-1](x), ans));
  plot([Phi, Theta, f], x=0..5, color=[red, blue, green]);

#
# No real domain assumption
#
  restart;
  eqn:= log[2](x^2 - 6*x) = 3 + log[2](1 - x);
  ans:=solve(eqn, x);
#
# x=2 gives the function complex values!
#
  simplify(eval(eqn, x=ans[1]));
  simplify(eval(eqn, x=ans[2]));

#
# Use real domain
#
  with(RealDomain):
  ans:=solve(eqn, x);
  simplify(eval(eqn, x=ans[1]));
  simplify(eval(eqn, x=ans[2]));