restart;
  with(plots):
  ic:=[seq(j, j=0..1, 0.1)]:
  display( [ seq
             ( odeplot
               ( dsolve
                 ( { diff(x(t),t)=y(t),
                     diff(y(t),t)=-x(t),
                     x(0)=1,
                     y(0)=ic[j]
                   },
                   numeric
                 ),
                 [[t, x(t)], [t, y(t)]],
                 t=0..10,
                 color=[red, blue]
               ),
               j=1..numelems(ic)
             )
           ]
         );

  restart;
  interface(rtablesize=20):
  with(StringTools):
  fname:="C:/Users/TomLeslie/Desktop/data.txt":
  T[0]:=Split( readline(fname), ";"):
  j:=1:
  line:=readline(fname):
  while line <> 0 do
        T[j]:=parse~( map2( Subs, ["E"="e", ","="."], Split(line, ";")));
        line:=readline(fname):
        j:=j+1;
  od:
  M:=Matrix([seq(T[k],k=0..j-1)]);

restart;
M:=Vector[column](10, [{{a, b}, {b, c}}, {{d, e}, {e, f}},
{{d, f}, {e, f}}, {{a, h}, {h, i}}, {{a, i}, {h, i}},
{{d, e}, {d, f}}, {{a, h}, {a, i}},{{c, d}}, {{a, b}, {a, c}}, {{a, c}, {b, c}}]):
Export("C:/Users/TomLeslie/Desktop/test.csv", M):

  restart;
#
# OP will have to change the file path in the following
# depnding on where (s)he, saved the "mpl" file
#
# Code in the file will read/execute with the following
# command, and produce a geometric construction, which
# happens to be a solution I provided for the problem
# posed here
#
# https://www.mapleprimes.com/questions/235455-Brianch0n-Theoreme-With-Hexagon
#
# The execution takes a few seconds, so wait
#
  read "C:/Users/TomLeslie/Desktop/hexplot3.mpl";

  restart;
  with(plots):
  local gamma:
  line := proc (x1, y1, x2, y2)
                options operator, arrow;
                (x-x1)*(y2-y1)-(y-y1)*(x2-x1)
          end proc:
  R := 3:
  ang := [0, (1/3)*Pi, 3*Pi*(1/4)+.2, 7*Pi*(1/6)+.4, 8*Pi*(1/5), 13*Pi*(1/7)]:
  for i to 6 do
      P || i := [R*cos(ang[i]), R*sin(ang[i])]
  end do:
  pts := [seq(P || i, i = 1 .. 6)]:
  for i to 6 do
      tang || i := x*P || i[1]+y*P || i[2] = R^2
  end do:
  sol:=solve({tang1,tang3},{x,y}):
  Q13:=[subs(sol,x),subs(sol,y)]:
  sol := solve({tang1, tang2}, {x, y}):
  Q1 := [subs(sol, x), subs(sol, y)]:
  sol := solve({tang2, tang3}, {x, y}):
  Q2 := [subs(sol, x), subs(sol, y)]:
  sol := solve({tang3, tang4}, {x, y}):
  Q3 := [subs(sol, x), subs(sol, y)]:
  sol := solve({tang4, tang5}, {x, y}):
  Q4 := [subs(sol, x), subs(sol, y)]:
  sol := solve({tang5, tang6}, {x, y}):
  Q5 := [subs(sol, x), subs(sol, y)]:
  sol := solve({tang1, tang6}, {x, y}):
  Q6 := [subs(sol, x), subs(sol, y)]:
  ptQ := [seq(Q || i, i = 1 .. 6)]:
  line14 := line(Q1[1], Q1[2], Q4[1], Q4[2]):
  L14 := implicitplot(line14, x = -10 .. 10, y = -10 .. 10, color = red):
  line25 := line(Q2[1], Q2[2], Q5[1], Q5[2]):
  L25 := implicitplot(line25, x = -10 .. 10, y = -10 .. 10, color = red):
  line36 := line(Q3[1], Q3[2], Q6[1], Q6[2]):
  L36 := implicitplot(line36, x = -10 .. 10, y = -10 .. 10, color = red):
  sol := solve({line14, line25}, {x, y}):
  I1 := [subs(sol, x), subs(sol, y)]:
  lineP23 := line(P2[1], P2[2], P3[1], P3[2]):
  lineP56 := line(P5[1], P5[2], P6[1], P6[2]):
  sol := solve({lineP23, lineP56}, {x, y}):
  gamma := [subs(sol, x), subs(sol, y)]:
  lineP12 := line(P1[1], P1[2], P2[1], P2[2]):
  lineP45 := line(P4[1], P4[2], P5[1], P5[2]):
  sol := solve({lineP12, lineP45}, {x, y}):
  beta := [subs(sol, x), subs(sol, y)]:
  lineP34 := line(P3[1], P3[2], P4[1], P4[2]):
  lineP16 := line(P1[1], P1[2], P6[1], P6[2]):
  sol := solve({lineP16, lineP34}, {x, y}):
  alpha := [subs(sol, x), subs(sol, y)]:
  pl:= line(alpha[1], alpha[2], gamma[1], gamma[2]):
#  hexa := seq(implicitplot(tang||i, x = -20 .. 20, y = -20 .. 20, linestyle=3,color = blue),i=1..6);
#  hexa:=plot([seq([P||i,P||(i mod 6)+1],i=1..6),color=green); 
  hexa:=plot([seq([`||`(P,i), `||`(P, irem(i,6)+1)],i=1..6)],thickness=4,color=green):
  hex:=plot([[Q1,Q2],[Q2,Q3],[Q3,Q4],[Q4,Q5],[Q5,Q6],[Q6,Q1]],thickness=4,color=green):
  tp := textplot({seq([op(pts[i]), cat("P", i)], i = 1 .. 6)}, 'align' = {'above', 'left'}):
  tpq := textplot({seq([op(ptQ[i]), cat("Q", i)], i = 1 .. 6)}, 'align' = {'above', 'left'}):
  TP:=textplot([[I1[],"I"],[alpha[],"alpha"],[beta[],"beta"],[gamma[],"gamma"]],'align' = {'above', 'left'}):
  slopes:=[seq(((dx,dy)->dy/dx)((pts[i]-pts[(i mod 6)+1])[]),i=1..6)]:
  lines:=zip((pt,slope)->y=slope*(x-pt[1])+pt[2],pts,slopes):
  plotpts:=plot(pts,style=point,colour=red,symbol=solidcircle,symbolsize=5):
  plotptQ:=plot(ptQ,style=point,colour=black,symbol=solidcircle,symbolsize=5):
  plotlines:=plot(rhs~(lines),style=line,linestyle=3,colour=magenta):
  cir:=implicitplot(x^2+y^2=R^2,x=-R..R,y=-R..R,color=black):
  PL:=implicitplot(pl,x=-1..15,y=-15..6,color=blue,thickness=3):
  display( plotpts,plotptQ,plotlines,hex,cir,L14,L25,L36,PL,tp,tpq,TP,
           axis = [gridlines = [4, color = blue]],
           view=[-6..10,-15..6],
           scaling=constrained,
           axes=none,
           size=[800,800]
         );

 

restart:
N[t]:=0.7: N[b]:=0.6:
Gr:=2.0: G[1]:=0.5: j[1]:=0.4: Pr:=2.5:
R:=1.3: Sc:=0.4:g[1]:=0.5:

F[0]:=0: F[1]:=1: F[2]:= A/2:
delta[-1]:=0: delta[0]:=1: delta[1]:=0:
Theta[0]:=1: Theta[1]:=B:
Phi[0]:=1: Phi[1]:=C:

for k from 2 to 20 do
    delta[k]:=0:
od:

for k from 0 to 3 do
    F[k + 3] := ( add((k - m + 1)*F[k - m + 1]*(m + 1)*F[m + 1], m = 0 .. k)*(1 + G[1]) + Gr*Theta[k] - add(F[k - m]*(m + 1)*(m + 2)*F[m + 2], m = 0 .. k) + (k + 1)*F[k + 1]*(M^2 + j[1]))/((k + 1)*(k + 2)*(k + 3)):
end do:

for k from 0 to 3 do
    Theta[k + 2] := -Pr*(N[t]*add(add(delta[l - 1]*(m - l + 1)*Theta[m - l + 1]*(k - m + 1)*Theta[k - m + 1], l = 0 .. m), m = 0 .. k) + add(F[k - m]*(m + 1)*Theta[m + 1], m = 0 .. k) + N[b]*add((k - m + 1)*Theta[k - m + 1]*(m + 1)*Phi[m + 1], m = 0 .. k) + Delta*Theta[k] - g[1]*add(Theta[k - m]*(m + 1)*F[m + 1], m = 0 .. k))/((1 + (4*R)/3)*(k + 1)*(k + 2)):
end do:

for k from 0 to 3 do
    Phi[k + 2] := -(Sc*add(F[k - m]*(m + 1)*Phi[m + 1], m = 0 .. k) + N[t]*(k + 1)*(k + 2)*Theta[k + 2]/N[b])/((k + 1)*(k + 2)):
end do:

f := add(F[i]*eta^i, i = 0 .. 3);
theta := add(Theta[i]*eta^i, i = 0 .. 3);
phi := add(Phi[i]*eta^i, i = 0 .. 3);

sol := RealDomain[solve]( { eval( diff(eval(f, [M = 1.1, Delta = 0.5]), eta), eta=3) = 1,
                            eval( theta, [M = 1.1, Delta = 0.5, eta=3]) = 0,
                            eval( phi, [M = 1.1, Delta = 0.5, eta=3]) = 0
                          }
                        )[3][];

d1 := diff( eval(f, [sol, M = 1.1, Delta = 0.5]), eta):
d2 := diff( eval(f, [sol, M = 1.6, Delta = 0.5]), eta):
d3 := diff( eval(f, [sol, M = 2.2, Delta = 0.5]), eta):
d4 := diff( eval(f, [sol, M = 2.6, Delta = 0.5]), eta):
plot([d1, d2, d3, d4], eta = 0 .. 7);

  s1 := eval(theta, [sol, Delta = 0.5] ):
  s2 := eval(theta, [sol, Delta = 1.0] ):
  s3 := eval(theta, [sol, Delta = 1.5] ):
  s4 := eval(theta, [sol, Delta = 2.0] ):
  plot([s1, s2, s3, s4], eta = 0 .. 7);