restart:

LumpSoliton_N4_TauFunction := proc()
    local i, j, eta, theta, B, f, exp34, exp3, exp4;

    # Step 1: Symbolic names
    for i from 1 to 4 do
        eta[i] := convert(cat("eta", i), name);
    end do:

    for i from 1 to 2 do
        theta[i] := convert(cat("theta", i), name);
    end do:

    for i from 1 to 3 do
        for j from i+1 to 4 do
            B[i,j] := convert(cat("B", i, j), name);
        end do:
    end do:

    # Step 2: Sub-expressions
    exp34 := exp(eta[3] + eta[4]) * (
        B[1,3]*B[2,3] + B[1,3]*B[2,4] + B[1,3]*theta[2]
      + B[1,4]*B[2,3] + B[1,4]*B[2,4] + B[1,4]*theta[2]
      + B[2,3]*theta[1] + B[2,4]*theta[1]
      + theta[1]*theta[2] + B[1,2]
    );

    exp3 := exp(eta[3]) * (
        B[1,3]*B[2,3] + B[1,3]*theta[2] + B[2,3]*theta[1] + theta[1]*theta[2] + B[1,2]
    );

    exp4 := exp(eta[4]) * (
        B[1,4]*B[2,4] + B[1,4]*theta[2] + B[2,4]*theta[1] + theta[1]*theta[2] + B[1,2]
    );

    # Step 3: Combine tau-function
    f := theta[1]*theta[2] + B[1,2] + exp3 + exp4 + A[3,4]*exp34;

    return simplify(f, size);
end proc:

LumpSoliton_N4_TauFunction();

restart:

LumpSolitonTauFunction := proc(N)
    local i, j, mu, mui, muval, eta, A, theta, B, f, mus, muvals, f0, expPart;

    # Create symbols eta[i], A[i,j]
    for i from 1 to N do
        eta[i] := convert(cat("eta", i), name);
    end do:
    for i from 1 to N-1 do
        for j from i+1 to N do
            A[i,j] := convert(cat("A", i, j), name);
        end do:
    end do:

    # === Case N = 1 ===
    if N = 1 then
        return exp(eta[1]);
    end if;

    # === Case N = 2: classic Hirota 2-soliton ===
    if N = 2 then
        return 1 + exp(eta[1]) + exp(eta[2]) + A[1,2]*exp(eta[1] + eta[2]);
    end if;

    # === Case N = 3: Lump + 1-soliton ===
    if N = 3 then
        for i from 1 to 2 do
            theta[i] := convert(cat("theta", i), name);
        end do:
        for i from 1 to 2 do
            for j from i+1 to 3 do
                B[i,j] := convert(cat("B", i, j), name);
            end do:
        end do:
        f := theta[1]*theta[2] + B[1,2]
           + exp(eta[3])*(B[1,3]*B[2,3] + theta[2]*B[1,3] + theta[1]*B[2,3] + theta[1]*theta[2] + B[1,2]);
        return simplify(f, size);
    end if;

    # === Case N = 4: Lump + 2-soliton (custom structure) ===
    if N = 4 then
        for i from 1 to 2 do
            theta[i] := convert(cat("theta", i), name);
        end do:
        for i from 1 to 3 do
            for j from i+1 to 4 do
                B[i,j] := convert(cat("B", i, j), name);
            end do:
        end do:
        f := theta[1]*theta[2] + B[1,2]
           + exp(eta[3])*(B[1,3]*B[2,3] + B[1,3]*theta[2] + B[2,3]*theta[1] + theta[1]*theta[2] + B[1,2])
           + exp(eta[4])*(B[1,4]*B[2,4] + B[1,4]*theta[2] + B[2,4]*theta[1] + theta[1]*theta[2] + B[1,2])
           + A[3,4]*exp(eta[3] + eta[4])*
                (B[1,3]*B[2,3] + B[1,3]*B[2,4] + B[1,3]*theta[2]
                + B[1,4]*B[2,3] + B[1,4]*B[2,4] + B[1,4]*theta[2]
                + B[2,3]*theta[1] + B[2,4]*theta[1] + theta[1]*theta[2] + B[1,2]);
        return simplify(f, size);
    end if;

    # === Case N > 4: General Hirota formula ===
    f0 := exp(add(mu[i]*eta[i], i = 1 .. N)
           + add(add(mu[i]*mu[j]*ln(A[i,j]), i = 1 .. j-1), j = 1 .. N)):
    mus := Iterator:-BinaryGrayCode(N):
    muvals := [seq({seq(mu[i] = mui[i], i = 1 .. N)}, mui in mus)]:
    f := add(eval(f0, muval), muval in muvals):

    # Clean up symbolic exp(...) terms
    with(StringTools):
    f := evalindets(f, specfunc(anything, exp),
        proc(expr)
            local symbols, Aterms, Etterms, Aprod, Et;
            symbols := indets(expr, name);
            Aterms, Etterms := selectremove(x -> evalb(Search("A", convert(x, string)) <> 0), symbols);
            Aprod := `if`(nops(Aterms) = 0, 1, mul(op(Aterms)));
            Et := add(op(Etterms));
            return Aprod * exp(Et);
        end proc
    );

    return simplify(f, size);
end proc:

LumpSolitonTauFunction(4);

LumpSolitonTauFunction(5);

restart:

LumpSolitonSuperposition_N3 := proc()
    local K, P, gamma, beta, theta, B, i, j, f, eta3;

    # Declare symbolic variables
    for i from 1 to 3 do
        K[i] := convert(cat("K", i), name);
        P[i] := convert(cat("P", i), name);
    end do:
    gamma := 'gamma': beta := 'beta':

    # Step 1: Define θ₁ and θ₂
    for i from 1 to 2 do
        theta[i] := simplify( - (K[i]^2 * gamma + P[i]^2 * beta)/K[i] + K[i]*y + P[i]*x );
    end do:

    # Step 2: Define B_{ij}
    B[1,2] := simplify( 6*P[1]*P[2]*K[1]*K[2]*(K[1]*P[1] + K[2]*P[2]) /
        ( (K[1]^2 - 2*P[1]) * (K[2]^2 - 2*P[2]) * K[1] * K[2] * (K[1]^2 * P[1]^2 + K[2]^2 * P[2]^2)) ):

    B[1,3] := simplify( 6*K[1]*P[1]*K[3]*P[3]*(K[1]*P[1] + K[3]*P[3]) /
        (-3*P[1]^2*P[3]*K[1]*K[3]^2*beta^2 - 2*P[1]*K[1]*K[3]^2*gamma + K[1]^2*P[1]^2) ):

    B[2,3] := simplify( 6*K[2]*P[2]*K[3]*P[3]*(K[2]*P[2] + K[3]*P[3]) /
        (-3*P[2]^2*P[3]*K[2]*K[3]^2*beta^2 - 2*P[2]*K[2]*K[3]^2*gamma + K[2]^2*P[2]^2) ):

    # Step 3: Define η₃ = eta3
    eta3 := convert("eta3", name):

    # Step 4: Build the tau-function f as in equation (13)
    f := theta[1]*theta[2] + B[1,2] + exp(eta3)*(B[1,3]*B[2,3] + theta[2]*B[1,3] + theta[1]*B[2,3] + theta[1]*theta[2] + B[1,2]):

    return simplify(f, size);
end proc:

LumpSolitonSuperposition_N3();

restart:

LumpSolitonTauFunction := proc(N)
    local i, j, mu, mui, muval, eta, A, f0, mus, muvals, f;
    local theta, B, expPart, lumpPart, eta3;

    # Define symbolic variables eta[i] and A[i,j]
    for i from 1 to N do
        eta[i] := convert(cat("eta", i), name);
    end do:
    for i from 1 to N-1 do
        for j from i+1 to N do
            A[i,j] := convert(cat("A", i, j), name);
        end do:
    end do:

    # Handle N = 1 separately (pure 1-soliton)
    if N = 1 then
        return exp(eta[1]);
    end if;

    # Handle N = 2: basic soliton interactions
    if N = 2 then
        return 1 + exp(eta[1]) + exp(eta[2]) + A[1,2]*exp(eta[1] + eta[2]);
    end if;

    # Handle N = 3: Lump + 1-soliton superposition
    if N = 3 then
        # Define θ₁ and θ₂
        for i from 1 to 2 do
            theta[i] := convert(cat("theta", i), name);
        end do:

        # Define B_{ij}
        for i from 1 to 2 do
            for j from i+1 to 3 do
                B[i,j] := convert(cat("B", i, j), name);
            end do:
        end do:

        eta3 := eta[3]:

        # Build symbolic expression like equation (13)
        lumpPart := theta[1]*theta[2] + B[1,2];
        expPart := exp(eta3) * (B[1,3]*B[2,3] + theta[2]*B[1,3] + theta[1]*B[2,3] + theta[1]*theta[2] + B[1,2]);
        f := lumpPart + expPart;

        return simplify(f, size);
    end if;

    # For N > 3, use generic Hirota expansion form
    f0 := exp(add(mu[i]*eta[i], i = 1 .. N)
           + add(add(mu[i]*mu[j]*ln(A[i,j]), i = 1 .. j-1), j = 1 .. N)):

    mus := Iterator:-BinaryGrayCode(N):
    muvals := [seq({seq(mu[i] = mui[i], i = 1 .. N)}, mui in mus)]:
    f := add(eval(f0, muval), muval in muvals):

    # Transform log terms into multiplicative form
    with(StringTools):
    f := evalindets(f, specfunc(anything, exp),
        proc(expr)
            local symbols, Aterms, Etterms, Aprod, Et;
            symbols := indets(expr, name);
            #Aterms, Etterms := selectremove(x -> Search("A", convert(x, string)) > 0, symbols);
             Aterms, Etterms := selectremove(x -> evalb(Search("A", convert(x, string)) <> 0), symbols);
            Aprod := if nops(Aterms) = 0 then 1 else mul(op(Aterms)) fi;
            Et := add(op(Etterms));
            return Aprod * exp(Et);
        end proc
    );

    return simplify(f, size);
end proc:

LumpSolitonTauFunction(1);

LumpSolitonTauFunction(2);

LumpSolitonTauFunction(3);

LumpSolitonTauFunction(4);

restart:

Get_f_N_eta_form := proc(N)
  local i, j, mu, eta, A, f0, mus, muvals, f,mui,muval ;

  # Step 1: Define eta[i] as symbolic variables (eta1, eta2, ...)
  for i from 1 to N do
    eta[i] := convert(cat('eta', i), name):
  end do:

  # Step 2: Define A[i,j] as symbolic variables (A12, A13, ...)
  for i from 1 to N-1 do
    for j from i+1 to N do
      A[i,j] := convert(cat('A', i, j), name):
    end do:
  end do:

  # Step 3: Define the generic Hirota exponential structure
  f0 := exp(add(mu[i]*eta[i], i = 1 .. N) +
            add(add(mu[i]*mu[j]*A[i,j], i = 1 .. j-1), j = 1 .. N)):

  # Step 4: Binary combinations of mu
  mus := Iterator:-BinaryGrayCode(N):
  muvals := seq({seq(mu[i] = mui[i], i = 1 .. N)}, mui in mus):

  # Step 5: Sum all evaluated cases
  f := add(eval(f0, muval), muval in muvals):

  return simplify(f);
end proc:

Get_f_N_eta_form(3);# output chanced  in another proc

restart:# manual

# Step 1: Number of solitons
N := 3:
printf("Step 1: N = %d\n", N);

# Step 2: Define eta[i] symbols
for i from 1 to N do
    eta[i] := convert(cat("eta", i), name);
end do:

# Step 3: Define A[i,j] symbols
for i from 1 to N-1 do
    for j from i+1 to N do
        A[i,j] := convert(cat("A", i, j), name);
    end do:
end do:

# Step 4: Build general f0 using log(Aij)
f0 := exp(add(mu[i]*eta[i], i = 1 .. N)
       + add(add(mu[i]*mu[j]*ln(A[i,j]), i = 1 .. j-1), j = 1 .. N)):

# Step 5: Generate all binary combinations
mus := Iterator:-BinaryGrayCode(N):
muvals := [seq({seq(mu[i] = mui[i], i = 1 .. N)}, mui in mus)]:

# Step 6: Evaluate sum over all combinations
f := add(eval(f0, muval), muval in muvals):

# Step 7: Load string tools and transform exp(A12 + eta1 + eta2) → A12 * exp(eta1 + eta2)
with(StringTools):

with(StringTools):

f := evalindets(f, specfunc(anything, exp),
    proc(expr)
        local symbols, Aterms, Etterms, Aprod, Et;
        symbols := indets(expr, name);
        
        # Separate Aij and eta terms properly
        Aterms, Etterms := selectremove(x -> Search("A", convert(x, string)) > 0, symbols);
        
        # Multiply Aij terms
        Aprod := if nops(Aterms) = 0 then 1 else mul(op(Aterms)) fi;
        
        # Sum eta terms
        Et := add(op(Etterms));
        
        return Aprod * exp(Et);
    end proc
);


# Step 8: Display cleaned result
printf("\nFinal tau function f_%d:\n", N);
f := simplify(f, size):
print(f);
 

Step 1: N = 3

Final tau function f_3:

restart:

TauFunction_Hirota := proc(N)
    local i, j, mu, mui, muval, eta, A, f0, mus, muvals, f;

    # Step 1: Define symbolic eta[i]
    for i from 1 to N do
        eta[i] := convert(cat("eta", i), name);
    end do:

    # Step 2: Define symbolic A[i,j]
    for i from 1 to N-1 do
        for j from i+1 to N do
            A[i,j] := convert(cat("A", i, j), name);
        end do:
    end do:

    # Step 3: Build the general exponential term f0
    f0 := exp(add(mu[i]*eta[i], i = 1 .. N)
           + add(add(mu[i]*mu[j]*ln(A[i,j]), i = 1 .. j-1), j = 1 .. N)):

    # Step 4: Generate all binary combinations of mu values
    mus := Iterator:-BinaryGrayCode(N):
    muvals := [seq({seq(mu[i] = mui[i], i = 1 .. N)}, mui in mus)]:

    # Step 5: Sum over all combinations
    f := add(eval(f0, muval), muval in muvals):

    # Step 6: Transform exp(Aij + eta terms) → Aij * exp(eta terms)
    f := evalindets(f, specfunc(anything, exp),
        proc(expr)
            local symbols, Aterms, Etterms, Aprod, Et;
            symbols := indets(expr, name);
            Aterms, Etterms := selectremove(x -> StringTools:-Search("A", convert(x, string)) > 0, symbols);
            Aprod := `if`(nops(Aterms) = 0, 1, mul(op(Aterms)));
            Et := add(op(Etterms));
            Aprod * exp(Et);
        end proc
    );

    # Step 7: Simplify and return the result
    f := simplify(f, size):
    return f;
end proc:

TauFunction_Hirota(3);

evalindets(f,`*`,
  proc(term)
    local inds:=indets(term,specindex(A));
    A[op(op~(inds))]*eval(term,inds=~1);
  end proc);

restart:

# Number of components
N := 2:

# Define ε-dependent variables
epsilon := 'epsilon':
for i from 1 to N do
    beta[i] := epsilon*theta[i];
end do:
B[1,2] := epsilon^2*Lambda12:

# Define the general exponent
f0 := exp(add(mu[i]*beta[i], i = 1 .. N) +
          add(add(mu[i]*mu[j]*B[i,j], i = 1 .. j-1), j = 1 .. N)):

# Generate all binary combinations of mu values
mus := Iterator:-BinaryGrayCode(N):
muvals := seq({seq(mu[i] = mui[i], i = 1 .. N)}, mui in mus):

# Sum over all binary cases
f := add(eval(f0, muval), muval in muvals):

# Expand as Taylor series in ε and extract ε² term (F2)
f_series := series(f, epsilon, 3):
F2 := simplify(coeff(f_series, epsilon, 2));  # This is F₂

restart:

GetF_N := proc(N)
  local i, j, epsilon, beta, theta, B, f0, f, f_series, mu, mus, muvals, F2, F2_clean, mui, muval;

  epsilon := 'epsilon':
  
  # Step 1: Define beta[i] = ε θ[i]
  for i from 1 to N do
    beta[i] := epsilon*theta[i]:
  end do:

  # Step 2: Define B[i,j] = ε² Λ[i,j]
  for i from 1 to N-1 do
    for j from i+1 to N do
      B[i,j] := epsilon^2*Lambda[i,j]:
    end do:
  end do:

  # Step 3: Build general Hirota term
  f0 := exp(add(mu[i]*beta[i], i = 1 .. N) +
            add(add(mu[i]*mu[j]*B[i,j], i = 1 .. j-1), j = 1 .. N)):

  # Step 4: Enumerate all binary mu combinations
  mus := Iterator:-BinaryGrayCode(N):
  muvals := seq({seq(mu[i] = mui[i], i = 1 .. N)}, mui in mus):

  # Step 5: Sum all evaluated exponential terms
  f := add(eval(f0, muval), muval in muvals):

  # Step 6: Expand in ε and extract ε² coefficient
  f_series := series(f, epsilon, 3):
  F2 := simplify(coeff(f_series, epsilon, 2)):

  # Step 7: Remove θ[i]^2 terms to get cross interactions only
  F2_clean := F2:
  for i from 1 to N do
    F2_clean := simplify(F2_clean - coeff(F2_clean, theta[i]^2)*theta[i]^2):
  end do:

  return F2_clean;
end proc:

GetF_N(2); # F

GetF_N(1);

GetF_N(3);

# version 1

restart:

# Load packages

with(plots):
local gamma;
# === Parameters === #
K1 := 'K1': P1 := 'P1':
K2 := 'K2': P2 := 'P2':
K3 := 'K3': P3 := 'P3':
beta := 'beta': gamma := 'gamma': alpha := 'alpha':
x := 'x': y := 'y': t := 't':

# === Theta functions (from equation (14)) === #
theta1 := (-K1^2*gamma + P1^2*beta)/K1 + K1*y + P1*x:
theta2 := (-K2^2*gamma + P2^2*beta)/K2 + K2*y + P2*x:

# === Soliton phase η³ === #
eta3 := K3*y + P3*x:

# === B_ij terms from equation (15) === #
B13 := beta * ((K1*K3*P1*(K1 + K3)) / (K1^2*P1^2 - P3^2) + K1*K3*P1*(K1 + K3)):
B23 := beta * ((K2*K3*P2*(K2 + K3)) / (K2^2*P2^2 - P3^2) + K2*K3*P2*(K2 + K3)):

# === Construct f (Equation 13) === #
f := theta1 + theta2 + exp(eta3)*(B13*B23 + theta2*B13 + theta1*B23 + theta1*theta2):

# === Final nonlinear wave field u === #
u := simplify(6/alpha * diff(ln(f), x$2)):

# === Optional: Plot u(x,y) for given parameters === #
# You can assign values and plot:
# K1 := 1: P1 := 2: etc...
# plot3d(eval(u, [K1=1, P1=2, ...]), x = -10 .. 10, y = -10 .. 10);

# version 2

restart:

# Load necessary packages
with(plots): with(VectorCalculus):
local  gamma;
# Declare symbolic parameters
K := [K1, K2, K3]: P := [P1, P2, P3]:
gamma := 'gamma': beta := 'beta': alpha := 'alpha':
x := 'x': y := 'y': t := 't':

# Define θ_i (Eq. 14): theta_i = (K_i^2 * γ + P_i^2 * β)/K_i + K_i*y + P_i*x
theta := [seq((K[i]^2*gamma + P[i]^2*beta)/K[i] + K[i]*y + P[i]*x, i = 1..2)]:

# Define η^3 = K3*y + P3*x + ω3*t + ξ3 (simplify to symbolic form)
eta3 := K[3]*y + P[3]*x + 'omega3'*t + 'xi3':

# Define B_ij for i<j<3 (Eq. 15)
B := Matrix(3, 3):
for i from 1 to 2 do
    for j from i+1 to 3 do
        B[i,j] := beta * (
            (6*K[i]*P[i]*K[j]*((K[i]*P[j] + K[j]*P[i])))/(
            -3*P[i]^2*K[i]^2 - 2*P[i]*K[i]*beta + beta*P[i]^2 + gamma*K[i]^2)
        );
    end do;
end do:

# Shortcuts for readability
theta1 := theta[1]: theta2 := theta[2]:
B12 := B[1,2]: B13 := B[1,3]: B23 := B[2,3]:

# Full expression of f from Eq. (13)
f := theta1 + theta2 + exp(eta3)*(B13*B23 + theta2*B13 + theta1*B23 + theta1*theta2 + B12):

# Compute u = (6/alpha) * ∂²/∂x² ln(f)
u := simplify(6/alpha * diff(ln(f), x$2)):

 

restart:

# Load packages
with(plots): with(VectorCalculus):

local gamma;
# Declare symbolic parameters
K := [K1, K2, K3, K4]: P := [P1, P2, P3, P4]:
gamma := 'gamma': beta := 'beta': alpha := 'alpha':
x := 'x': y := 'y': t := 't':

# Define θ_i terms (Eq. 14)
theta := [seq((K[i]^2*gamma + P[i]^2*beta)/K[i] + K[i]*y + P[i]*x, i = 1..2)]:
theta1 := theta[1]: theta2 := theta[2]:

# Define η^3 and η^4
eta3 := K[3]*y + P[3]*x + 'omega3'*t + 'xi3':
eta4 := K[4]*y + P[4]*x + 'omega4'*t + 'xi4':

# Define B_ij matrix using general formula from Eq. (15)
B := Matrix(4, 4):
for i from 1 to 4 do
    for j from i+1 to 4 do
        B[i,j] := beta * (
            (6*K[i]*P[i]*K[j]*((K[i]*P[j] + K[j]*P[i])))/(
            -3*P[i]^2*K[i]^2 - 2*P[i]*K[i]*beta + beta*P[i]^2 + gamma*K[i]^2)
        );
    end do;
end do:

# Shortcuts for B_ij terms
B12 := B[1,2]: B13 := B[1,3]: B23 := B[2,3]:
B14 := B[1,4]: B24 := B[2,4]:

# Assume symbolic A34 (interaction coefficient)
A34 := 'A34':

# Construct f as in Equation (27)
f := theta1 + theta2 + B12
     + A34 * exp(eta3 + eta4) * (B13*B23 + B13*theta2 + B23*theta1 + theta1*theta2 + B12)
     + B13*theta2 + B12*B23 + B14*B24 + B14*theta2 + B24*theta1 + B24*theta1*theta2
     + (B13*B23 + B12) * exp(eta3)
     + (B14*B24 + B12) * exp(eta4):

# Compute u = (6/alpha) * ∂²/∂x² ln(f)
u := simplify(6/alpha * diff(ln(f), x$2)):


 

restart:

with(plots): with(VectorCalculus):

# === Parameters === #
K := [K1, K2, K3, K4, K5]: P := [P1, P2, P3, P4, P5]:
gamma := 'gamma': beta := 'beta': alpha := 'alpha':
x := 'x': y := 'y': t := 't':

# === θ_i terms for lump wave === #
theta := [seq((K[i]^2*gamma + P[i]^2*beta)/K[i] + K[i]*y + P[i]*x, i = 1..2)]:
theta1 := theta[1]: theta2 := theta[2]:

# === η^j soliton phases === #
eta := [seq(K[i]*y + P[i]*x + 'omega'||i*t + 'xi'||i, i=3..5)]:

# === Define B_ij terms === #

B := Matrix(5, 5):
for i from 1 to 5 do
    for j from i+1 to 5 do
        B[i, j] := beta * (
            (6*K[i]*P[i]*K[j]*(K[i]*P[j] + K[j]*P[i])) /
            (-3*P[i]^2*K[i]^2 - 2*P[i]*K[i]*beta + beta*P[i]^2 + gamma*K[i]^2)
        );
    end do;
end do:

# === A_ij interaction terms — symbolic form === #
A34 := 'A34': A35 := 'A35': A45 := 'A45':

# === Reuse B_ij notation === #
B12 := B[1,2]: B13 := B[1,3]: B14 := B[1,4]: B15 := B[1,5]:
B23 := B[2,3]: B24 := B[2,4]: B25 := B[2,5]:
B34 := B[3,4]: B35 := B[3,5]: B45 := B[4,5]:

# === Build f (Eq. 34) === #
f := theta1 + theta2 + B12
  + (B13 + theta1)*(B23 + theta2)*exp(eta[1])
  + (B14 + theta1)*(B24 + theta2)*exp(eta[2])
  + (B15 + theta1)*(B25 + theta2)*exp(eta[3])
  + A34*(B13 + theta1)*(B23 + theta2)*exp(eta[1] + eta[2])
  + A35*(B13 + theta1)*(B23 + theta2)*exp(eta[1] + eta[3])
  + A45*(B14 + theta1)*(B24 + theta2)*exp(eta[2] + eta[3])
  + A34*A35*A45*(B13 + B14 + B15)*(B23 + B24 + B25)*exp(eta[1] + eta[2] + eta[3]):

# === Final nonlinear wave field u === #
u := simplify(6/alpha * diff(ln(f), x$2)):