restart:
with(plots):

### Lump Parameters ### PROCEDURE ZONDER INITIELE POSITIE LUMPS ?
a := -1; b := -2; alpha := 1; beta := 1;
vx := -3/5; vy := 1/5;
trajectory_slope := vy/vx;  # automatisch bepaald
U := (x, y, t, a, b, alpha, beta) ->
    4*beta*((a*y - 2*alpha*t + x)*b^2 + a*(-2*beta*t + a*(a*y - 2*alpha*t + x))*b^2) /
    (-b^6*beta*y^2
    + (-4*t*y*beta^2 + (-2*a^2*y^2 + (4*alpha*t - 2*x)*a*y - 4*(alpha*t - x/2)^2)*beta + 3*a)*b^4
    + (-4*t^2*beta^3 + 4*a*t*(a*y - 2*alpha*t + x)*beta^2 - a^2*(a*y - 2*alpha*t + x)^2*beta + 6*a^3)*b^2
    + 3*a^5);

### �� Dynamische Lump-Procedure (met parameter-uitleg) �� ###
DynamicLumpPlot := proc(input::list, mode::symbol := 'tijd')
  local i, n, pos_initial, tijden, pos_shift, contour, kleuren, trajectory_plot, text_labels, U, extra_info, xmin, xmax, ymin, ymax, param_info;

  n := nops(input);
  pos_initial := [seq([200*(i-(n+1)/2), -200*(i-(n+1)/2)/3], i=1..n)];

  if mode = 'tijd' then
    tijden := input;
    for i from 1 to n do
      pos_shift[i] := [pos_initial[i][1] + vx*tijden[i], pos_initial[i][2] + vy*tijden[i]];
    end do:
    extra_info := cat(seq(sprintf("t=%.2f: (%.2f, %.2f)   ", tijden[i], pos_shift[i][1], pos_shift[i][2]), i=1..n));

  elif mode = 'positie' then
    pos_shift := input;
    for i from 1 to n do
      tijden[i] := (pos_shift[i][1] - pos_initial[i][1]) / vx;
      if abs(pos_shift[i][2] - trajectory_slope*pos_shift[i][1]) > 1e-3 then
        printf("⚠️ Lump %d ligt NIET precies op trajectory (afwijking=%.2f). Correcte positie: (%.2f, %.2f)\n",
          i, abs(pos_shift[i][2] - trajectory_slope*pos_shift[i][1]), pos_shift[i][1], trajectory_slope*pos_shift[i][1]);
      end if;
    end do:
    extra_info := cat(seq(sprintf("(%.2f, %.2f): t=%.2f   ", pos_shift[i][1], pos_shift[i][2], tijden[i]), i=1..n));

  else
    error "Kies modus 'tijd' of 'positie'";
  end if:

  xmin := min(seq(pos_shift[i][1], i=1..n)) - 50; xmax := max(seq(pos_shift[i][1], i=1..n)) + 50;
  ymin := min(seq(pos_shift[i][2], i=1..n)) - 50; ymax := max(seq(pos_shift[i][2], i=1..n)) + 50;

  xmin := min(xmin, -200): xmax := max(xmax, 200):
  ymin := min(ymin, trajectory_slope*xmin - 50): ymax := max(ymax, trajectory_slope*xmax + 50):

  U := (x, y, t, a, b, alpha, beta) ->
    4*beta*((a*y - 2*alpha*t + x)*b^2 + a*(-2*beta*t + a*(a*y - 2*alpha*t + x))*b^2) /
    (-b^6*beta*y^2
    + (-4*t*y*beta^2 + (-2*a^2*y^2 + (4*alpha*t - 2*x)*a*y - 4*(alpha*t - x/2)^2)*beta + 3*a)*b^4
    + (-4*t^2*beta^3 + 4*a*t*(a*y - 2*alpha*t + x)*beta^2 - a^2*(a*y - 2*alpha*t + x)^2*beta + 6*a^3)*b^2
    + 3*a^5);

  kleuren := [red, blue, green, orange, purple, cyan, magenta]:

  for i from 1 to n do
    contour[i] := contourplot(
      eval(U(x - pos_shift[i][1], y - pos_shift[i][2], 0, a, b, alpha, beta)),
      x = xmin..xmax, y = ymin..ymax, contours=30, color=kleuren[i mod 7 + 1],
      grid=[200,200], transparency=0.1
    );
  end do:

  trajectory_plot := implicitplot(y=trajectory_slope*x, x=xmin..xmax, y=ymin..ymax, color=black, thickness=2):

  text_labels := textplot([
    seq([pos_shift[i][1], pos_shift[i][2]+(ymax-ymin)*0.03,
         cat("t=", evalf[3](tijden[i])),
         'font'=["Arial",12,'bold'], 'color'=black], i=1..n)
  ]):

  param_info := sprintf("Lump-parameters betekenis:\n a=%.2f (richting/vormfactor), b=%.2f (breedte/hoogte-factor), α=%.2f, β=%.2f (amplitude/richting)\n → Trajectory-slope = %.2f",
                         a,b,alpha,beta,trajectory_slope);

  display(
    [seq(contour[i], i=1..n)], trajectory_plot, text_labels,
    title = cat(n, "-Lump Soliton(en) Dynamisch"),
    labels = ["x","y"], scaling = constrained,
    caption = typeset(extra_info, "\n\n", param_info)
  );
end proc:

DynamicLumpPlot([0,1,2,3,4,5,6,7,8,9,10], 'tijd');# time in sec ?

DynamicLumpPlot([[-150,50],[-80,80/3],[50,-50/3],[0,0]],'positie');# position

DynamicLumpPlot([[-150,50]],'positie');