restart;

# Definities
g := a1*x + a2*y + a3*z + a4*t + a5:
h := a6*x + a7*y + a8*z + a9*t + a10:
f := g^2 + h^2 + a11:

# Afgeleiden
fx := diff(f, x): fy := diff(f, y): fz := diff(f, z): ft := diff(f, t):
fxx := diff(fx, x): fxy := diff(fx, y): fxz := diff(fx, z): fyy := diff(fy, y):
fxxx := diff(fxx, x): fxxy := diff(fxx, y): fxxxy := diff(fxxx, y):
fxxxx := diff(fxxx, x): fxxxxx := diff(fxxxx, x): fxxxxxx := diff(fxxxxx, x):

# PDE
PDE := 9*(f*diff(f, x, t) - fx*ft)
    -5*(f*diff(f, x, x, x, y) - 3*fx*diff(f, x, x, y) + 3*fxx*fxy - fxxx*fy)
    + (f*diff(f, x, x, x, x, x, x) - 6*fx*diff(f, x, x, x, x, x) + 15*fxx*diff(f, x, x, x, x) - 10*(fxxx)^2)
    -5*(f*fyy - fy^2)
    +alpha*(f*fxx - fx^2)
    +beta*(f*fxy - fx*fy)
    +gamma*(f*fxz - fx*fz):

PDE := expand(PDE): PDE := simplify(PDE):

# Coëfficiënten pakken
coeff1 := coeff(PDE, x, 1):
coeff2 := coeff(PDE, y, 1):
coeff3 := coeff(PDE, z, 1):
coeff4 := coeff(PDE, t, 1):
coeff_const := eval(PDE, {x=0, y=0, z=0, t=0}):

# Stelsel van vergelijkingen
Eqns := {coeff1=0, coeff2=0, coeff3=0, coeff4=0, coeff_const=0}:
# Stelsel
##Eqns := {coeff_x=0, coeff_y=0, coeff_z=0, coeff_t=0, coeff_const=0}:

# Oplossen
Sol := solve(Eqns, {a4, a9, a11}):
Sol;

restart:

DrawOrderedLines := proc(shape1::list, isolated1::list, shape2::list, isolated2::list)
  local i, p1, p2, allplots, colorlist, polyline, isolateplot, j:
  uses plots, plottools;
  colorlist := [blue, red]:
  allplots := []:

  for i from 1 to 2 do
    if i = 1 then
      polyline := shape1:
      isolateplot := isolated1:
    else
      polyline := shape2:
      isolateplot := isolated2:
    end if:

    # Verbind punten in exacte volgorde
    if nops(polyline) >= 2 then
      for j from 1 to nops(polyline) - 1 do
        p1 := polyline[j]:
        p2 := polyline[j+1]:
        allplots := [op(allplots), line(p1, p2, color=colorlist[i], thickness=2)]:
      end do:
    end if:

    # Teken losse geïsoleerde punten
    if nops(isolateplot) > 0 then
      allplots := [op(allplots), pointplot(isolateplot, color=colorlist[i], symbol=solidcircle, symbolsize=15)]:
    end if:
  end do:

  return display(allplots, scaling=constrained, title="");
end proc:

fig1 := [[0, 0], [4, 0], [4, 3], [0, 3], [5, 5],[8,1],[8,8],[-5,8]]:
isolated1 := [[0,3],[2, 1],[-5,8]]:
fig2 := [[6, 6], [8, 6], [8, 8], [6, 8], [9, 9]]:
isolated2 := [[7, 7]]:

DrawOrderedLines(fig1, isolated1, fig2, isolated2);

restart:
with(DEtools):

ExactODEGroupingInteractive := proc(ODE, x, y)
  local lhs_expr, M, N, M_orig, N_orig, dM_dy, dN_dx, matched,
        intM, u_partial, du_dy, gprime, g, u_final, terms,
        dx_term, dy_term, ODE_classical, advice;

  print("�� STEP 0: Start by reading the differential equation in differential form:");
  lhs_expr := lhs(ODE);
  print("   → Equation provided: ", ODE);

  print("\n�� STEP 0B: Try to classify the ODE using Maple’s odeadvisor");

  # Extract M and N
  terms := [op(lhs_expr)];
  dx_term := select(t -> has(t, Dx), terms)[1];
  dy_term := select(t -> has(t, Dy), terms)[1];
  M := simplify(dx_term / Dx);
  N := simplify(dy_term / Dy);

  # Save for later
  M_orig := M:
  N_orig := N:

  # Substitute y = y(x)
  M := subs(y = y(x), M):
  N := subs(y = y(x), N):

  print("   We now rewrite the ODE in the form:");
  print("     M(x,y)*dy/dx + N(x,y) = 0");
  ODE_classical := M + N*diff(y(x), x) = 0;
  print("   → Rewritten ODE: ", ODE_classical);

  print("   Ask Maple’s odeadvisor for classification:");
  advice := odeadvisor(ODE_classical, y(x));
  print("�� Maple's odeadvisor says: ", advice);

  if member(_exact, advice) then
    print("✅ Great! Maple recognizes this ODE as exact.");
  else
    print("⚠️ Maple did NOT classify this as exact.");
    print("   → We will verify exactness ourselves using partial derivatives.");
  end if;

  print("\n�� STEP 1: Extract the functions M(x, y) and N(x, y) from the equation:");
  print("   → M(x, y) = ", M_orig);
  print("   → N(x, y) = ", N_orig);

  print("\n�� STEP 2: Check if the ODE is exact by comparing partial derivatives:");
  print("   → Compute ∂M/∂y and ∂N/∂x:");
  dM_dy := diff(M_orig, y):
  dN_dx := diff(N_orig, x):
  print("     ∂M/∂y = ", dM_dy);
  print("     ∂N/∂x = ", dN_dx);
  matched := simplify(dM_dy - dN_dx);
  if matched <> 0 then
    return "❌ This equation is NOT exact. Grouping method is not applicable.";
  else
    print("✅ Yes! The equation is exact.");
  end if;

  print("\n�� STEP 3: Integrate M(x, y) with respect to x (treat y as constant):");
  print("   → ∫ M(x, y) dx = ");
  intM := int(M_orig, x):
  print("     = ", intM, " + g(y)");

  print("\n�� STEP 4: Differentiate that result with respect to y:");
  print("   → ∂/∂y of the integral:");
  du_dy := diff(intM, y):
  print("     ∂(partial result)/∂y = ", du_dy);

  print("\n�� STEP 5: Now solve for g'(y):");
  print("   → g'(y) = N(x, y) - ∂/∂y(previous result)");
  gprime := simplify(N_orig - du_dy):
  print("     g'(y) = ", gprime);

  print("\n�� STEP 6: Integrate g'(y) with respect to y:");
  print("   → ∫ g'(y) dy = ");
  g := int(gprime, y):
  print("     g(y) = ", g);

  print("\n�� STEP 7: Add both parts to find the general solution u(x, y):");
  u_final := simplify(intM + g):
  print("     u(x, y) = ", u_final);

  print("\n�� FINAL RESULT: The general solution to the ODE is:");
  return u_final = 'C';

end proc:

# Example call
ExactODEGroupingInteractive(
  (6*x*y^2 + 4*x^3*y)*Dx + (6*x^2*y + x^4 + exp(y))*Dy = 0,
  x, y);

ExactODEGroupingInteractive(
  (3*x*y^2 + 4*x^3*y)*Dx + (6*x^2*y + x^4 + exp(y))*Dy = 0,
  x, y);

complex_map_plot := proc(f, zcurveLijst, t_range,
                         {kleurenset := ["blue","red","green","magenta","orange","cyan"],
                          lijnstijl := solid, resolutie := 200,
                          domein_naam := "Krommes in z-vlak", toon_label := false})
    local i, t, zcurve, ReZ, ImZ, ReW, ImW, zplots, wplots, kleur, zplot, wplot, annotatie, labeltekst;
    uses plots, plottools, InertForm;

    zplots := [];  wplots := [];

    for i to nops(zcurveLijst) do
        zcurve := zcurveLijst[i];
        kleur := kleurenset[i mod nops(kleurenset) + 1];

        ReZ := t -> Re(zcurve(t));
        ImZ := t -> Im(zcurve(t));
        ReW := t -> Re(f(zcurve(t)));
        ImW := t -> Im(f(zcurve(t)));

        zplot := plot([ReZ(t), ImZ(t), t = t_range],
                      color = kleur, linestyle = lijnstijl, thickness = 2,
                      scaling = constrained, axes = boxed, numpoints = resolutie,
                      labels = ["Re(z)", "Im(z)"]);

        if toon_label then
            labeltekst := typeset("z(t) = ", InertForm:-Display(zcurve(t)));
            annotatie := textplot([ReZ((rhs(t_range) + lhs(t_range))/2),
                                   ImZ((rhs(t_range) + lhs(t_range))/2),
                                   labeltekst], font = [TIMES, BOLD, 12]);
            zplot := display(zplot, annotatie);
        end if;

        wplot := plot([ReW(t), ImW(t), t = t_range],
                      color = kleur, linestyle = lijnstijl, thickness = 2,
                      scaling = constrained, axes = boxed, numpoints = resolutie,
                      labels = ["Re(w)", "Im(w)"]);

        zplots := [op(zplots), zplot];
        wplots := [op(wplots), wplot];
    end do;

    display(Array([display(zplots, title = cat("z-vlak: ", domein_naam)),
                   display(wplots, title = "w-vlak: Beeld onder f(z)")]), insequence = false);
end proc:

f := z -> exp(z);
zcurves := [t -> -1 + I*t, t -> 0 + I*t, t -> 1 + I*t, t -> 2 + I*t];
complex_map_plot(f, zcurves, -Pi..Pi, toon_label = false, domein_naam = "Verticale lijnen Re(z) = -1,0,1,2");

f := z -> exp(z);
z_driehoek := proc(t) piecewise(
  t <= 1, 2*t,
  t <= 2, 2 + (t-1)*(-1+I),
  t <= 3, 1 + I + (t-2)*(-1 - I)
) end proc;

complex_map_plot(f, [z_driehoek], 0..3, toon_label=true, domein_naam="Driehoek");

zcurves := [
    t -> t - 2*I,
    t -> t - I,
    t -> t,
    t -> t + I,
    t -> t + 2*I
];

f := z -> exp(z);
complex_map_plot(f, zcurves, -2..2, domein_naam = "Horizontale lijnen Im(z) = -2..2", toon_label = false);

restart;
with(plots):
with(DEtools):

SEIRTPplot := proc(param::list, tmin::numeric, tmax::numeric, dynRange::boolean := true, yRange::range := 0..1200)
    description "SEIRTP Model Plotter – Simulates and visualizes the effects of parameter changes on all compartments";

    local Lambda, alpha, delta, K, r, tau, gamma_, mu, phi, eta, sigma, beta, theta, xi;
    local S, E, Infectious, R, P, T, t, de, init, sol, plt, p;

    # Default parameter values
    Lambda := 0.0133:   alpha := 0.07954551: delta := 0.9: K := 300:
    r := 0.076:         tau := 0.0900982:   gamma_ := 0.00917: mu := 0.00056:
    phi := 0.09:        eta := 0.009:       sigma := 0.000456: beta := 0.567:
    theta := 0.009:     xi := 0.0487:

    # Overwrite defaults if needed
    for p in param do eval(p); end do;

    # Differential equations
    de := [
      diff(S(t),t) = -P(t)*S(t)*alpha - S(t)*mu + Lambda,
      diff(E(t),t) = alpha*S(t)*P(t) - (-T(t)*eta + 1)*beta*E(t) - theta*E(t) - mu*E(t),
      diff(Infectious(t),t) = (-T(t)*eta + 1)*beta*E(t) - delta*Infectious(t) - gamma_*Infectious(t) - mu*Infectious(t),
      diff(R(t),t) = theta*E(t) + gamma_*Infectious(t) - mu*R(t),
      diff(P(t),t) = xi*Infectious(t) - sigma*T(t)*P(t) - tau*P(t),
      diff(T(t),t) = r*T(t)*(1 - T(t)/K) - phi*T(t)
    ]:

    # Initial conditions
    init := [S(0)=1000, E(0)=10, Infectious(0)=5, R(0)=0, P(0)=1, T(0)=1]:

    # Numerical solution
    sol := dsolve({op(de), op(init)}, numeric, output = listprocedure, range=tmin..tmax):

    # Plot generation
    if dynRange then
        plt := odeplot(sol,
            [[t, S(t)], [t, E(t)], [t, Infectious(t)], [t, R(t)], [t, P(t)], [t, T(t)]],
            tmin..tmax,
            legend = ["Susceptible S(t)", "Exposed E(t)", "Infectious I(t)", "Recovered R(t)", "Pathogen P(t)", "Treatment T(t)"],
            title = "SEIRTP Model – Dynamic Axis Scaling",
            thickness = 2,
            style = [line]);
    else
        plt := odeplot(sol,
            [[t, S(t)], [t, E(t)], [t, Infectious(t)], [t, R(t)], [t, P(t)], [t, T(t)]],
            tmin..tmax,
            legend = ["Susceptible S(t)", "Exposed E(t)", "Infectious I(t)", "Recovered R(t)", "Pathogen P(t)", "Treatment T(t)"],
            title = cat("SEIRTP Model – Fixed Axis Scaling ", convert(yRange,string)),
            thickness = 2,
            style = [line],
            view = [tmin..tmax, yRange]);
    end if;

    # �� Procedure Input Explanation
    printf~([
      "��️ Procedure Input Explanation:\n",
      sprintf("%-10s : %s\n", "param",    "List of parameter assignments (e.g. [r = 0.12, beta = 0.6])"),
      sprintf("%-10s : %s\n", "tmin",     "Start time for simulation (e.g. 0)"),
      sprintf("%-10s : %s\n", "tmax",     "End time for simulation (e.g. 100)"),
      sprintf("%-10s : %s\n", "dynRange", "true = dynamic y-axis scaling, false = use yRange"),
      sprintf("%-10s : %s\n", "yRange",   "Optional fixed y-axis range (e.g. 0..1500)")
    ]);

    # �� Separation line
    printf("------------------------------------------------------------\n");

    # �� Parameter Meaning
    printf~([
      "�� Parameter Meaning:\n",
      sprintf("%-10s : %s\n", "Lambda",   "Recruitment rate into S"),
      sprintf("%-10s : %s\n", "alpha",    "Infection rate via environment (S to E)"),
      sprintf("%-10s : %s\n", "beta",     "Progression rate from E to I"),
      sprintf("%-10s : %s\n", "delta",    "Disease-induced death rate"),
      sprintf("%-10s : %s\n", "gamma_",   "Recovery rate from I"),
      sprintf("%-10s : %s\n", "mu",       "Natural death rate"),
      sprintf("%-10s : %s\n", "theta",    "Direct recovery rate from E"),
      sprintf("%-10s : %s\n", "xi",       "Pathogen production by I"),
      sprintf("%-10s : %s\n", "sigma",    "Treatment effectiveness against P"),
      sprintf("%-10s : %s\n", "tau",      "Natural pathogen decay"),
      sprintf("%-10s : %s\n", "r",        "Growth rate of treatment capacity"),
      sprintf("%-10s : %s\n", "K",        "Carrying capacity of treatment system"),
      sprintf("%-10s : %s\n", "phi",      "Treatment dropout rate"),
      sprintf("%-10s : %s\n", "eta",      "Effect of treatment on infection process")
    ]);

    return plt;
end proc:

# Dynamisch bereik
SEIRTPplot([], 0, 15, true);
 

��️ Procedure Input Explanation:
param      : List of parameter assignments (e.g. [r = 0.12, beta = 0.6])
tmin       : Start time for simulation (e.g. 0)
tmax       : End time for simulation (e.g. 100)
dynRange   : true = dynamic y-axis scaling, false = use yRange
yRange     : Optional fixed y-axis range (e.g. 0..1500)
------------------------------------------------------------
�� Parameter Meaning:
Lambda     : Recruitment rate into S
alpha      : Infection rate via environment (S to E)
beta       : Progression rate from E to I
delta      : Disease-induced death rate
gamma_     : Recovery rate from I
mu         : Natural death rate
theta      : Direct recovery rate from E
xi         : Pathogen production by I
sigma      : Treatment effectiveness against P
tau        : Natural pathogen decay
r          : Growth rate of treatment capacity
K          : Carrying capacity of treatment system
phi        : Treatment dropout rate
eta        : Effect of treatment on infection process

# Vast bereik met aangepaste schaal
SEIRTPplot([r = 0.12], 0, 15, false, 0..1200);

��️ Procedure Input Explanation:
param      : List of parameter assignments (e.g. [r = 0.12, beta = 0.6])
tmin       : Start time for simulation (e.g. 0)
tmax       : End time for simulation (e.g. 100)
dynRange   : true = dynamic y-axis scaling, false = use yRange
yRange     : Optional fixed y-axis range (e.g. 0..1500)
------------------------------------------------------------
�� Parameter Meaning:
Lambda     : Recruitment rate into S
alpha      : Infection rate via environment (S to E)
beta       : Progression rate from E to I
delta      : Disease-induced death rate
gamma_     : Recovery rate from I
mu         : Natural death rate
theta      : Direct recovery rate from E
xi         : Pathogen production by I
sigma      : Treatment effectiveness against P
tau        : Natural pathogen decay
r          : Growth rate of treatment capacity
K          : Carrying capacity of treatment system
phi        : Treatment dropout rate
eta        : Effect of treatment on infection process