differential equations – NDSolve not recognizing symbolic derivative as a function

I’m trying to solve a nonlinear first order initial value problem ${x'(t)= f(x,t), x(t_0)=x_0 }$ using NDSolve. The function $f(x,t)$ is complicated, as is the computation of the initial x-value. The code computes these values accurately, but when I go to solve the IVP I get the error message:

NDSolve::ndode: The equations {True,0==SolSpeed(xsol(tsol)/tsol,0.25,1,2)} are not differential equations or initial conditions in the dependent variables {xsol}.

I’ve read lots of help files and tried restarting the kernel, but nothing seems to help. I’m sure my code isn’t optimal, but I can’t figure out why it won’t run at all.

Any ideas?

Here is the full code, the NDSolve call is at the bottom inside a function called SolPeakLeading.

XiCrit(eta1_, eta2_) :=
 Module({Whitham, m},
  Whitham(m_?NumberQ) = 
   4 (1 - m) EllipticK(m)/EllipticE(m) + 2 (1 + m);
  eta2^2*Whitham(eta1^2/eta2^2)
  )

ComputeAlpha(xi_?NumericQ, eta1_?NumericQ, eta2_?NumericQ) :=
 
 Module({xi$critical},
  xi$critical = XiCrit(eta1, eta2);
  If(xi >= xi$critical,
   eta2,
   Module({Whitham, asq, ans, m},
    Whitham(m_?NumberQ) = 
     4 (1 - m) EllipticK(m)/EllipticE(m) + 2 (1 + m);
    asq = 
     asq /. FindRoot(
       xi == asq *
         Whitham(eta1^2/asq), {asq, ( 
          eta1 + (eta2 - 
              eta1)*(xi - 4 eta1^2)/(xi$critical - 4 eta1^2))^2},
       MaxIterations -> 150
       );
    ans = Sqrt(asq);
    (* Now test if FindRoot gives a spurious complex answer *)
    
    If(Im (ans) != 0,
     (* If a nonzero real part switch to a bracketing method to 
compute alpha *)
     Clear(asq);
     Sqrt(
      asq /. FindRoot(
        xi == asq*Whitham(eta1^2/asq), {asq, eta1^2, eta2^2}, 
        Method -> "Secant")),
     (* If false, return the original answer *)
     ans
     )
    )
   )
  )

sol(xval_, tval_, eta1_, eta2_, kappa_, x0_, (Sigma)_) :=
 
 If(xval < 4 eta1^2*tval,
  2 (Sigma) kappa Sech(2 kappa (xval - x0 - 4 kappa^2 tval)), (* 
  If true, evaluate the zero background soliton *)
  (* If false, 
  the module computes the background and soliton solution in the 
elliptic regions *)
  Module(
   {
    xi, alpha, R, m1, eK, eK1, nome, (Gamma)sq, Abel, (CurlyPhi), 
    snSlow, dnSlow, cnSlow,
    Xi, TT, Omega, qBackground, Qminus, QminusPrime, w22sq, Y, X, 
    qsol
    },
   xi = xval/tval;
   alpha = ComputeAlpha(xi, eta1, eta2); 
   R = -Sqrt(( kappa^2 - alpha^2) (kappa^2 - eta1^2)); (* 
   R(k) evaluated at k=i(Kappa) *)
   
   m1 = 4*alpha*eta1/(alpha + eta1)^2;
   eK = EllipticK(eta1^2/alpha^2);
   eK1 = EllipticK(m1);
   nome = Exp(-Pi EllipticK(1 - eta1^2/alpha^2)/(2 eK)) ;
   (Gamma)sq = 
    Sqrt((kappa - alpha)/(kappa - eta1) * (kappa + eta1)/(kappa + 
         alpha)); (* (Gamma)^2 evaluated at k=i(Kappa) *)
   
   Abel = -InverseJacobiDC(kappa/alpha, eta1^2/alpha^2)/(4 eK); (* 
   Abel evaluated at k=i(Kappa) > i (Alpha)*)
   
   snSlow = JacobiSN(2 eK1 (Abel + 1/4), m1);
   cnSlow = JacobiCN(2 eK1 (Abel + 1/4), m1);
   dnSlow = JacobiDN(2 eK1 (Abel + 1/4), m1);
   (* 
   Everything before this is a function of the slow evolution 
variable (Xi). 
   Everything after depends on the fast (order one) evolution in xval 
and tval
   *)
   Omega(Xi_, TT_) = -Pi*alpha/eK*TT*(Xi - 2 (eta1^2 + alpha^2));
   (* Compared with notes (CurlyPhi) (here) is really -
   I*(CurlyPhi) + log((Chi))/
   2 (from notes)  *)
   (CurlyPhi)(Xi_, TT_) = 
    R *TT*(4 kappa - (1/
           kappa)*(EllipticPi(eta1^2/kappa^2, eta1^2/alpha^2)/
           eK)*(Xi - 2 (eta1^2 + alpha^2))) - kappa x0 + 
     Log(2 kappa)/2;
   qBackground(Xi_, 
     TT_) = (alpha + 
       eta1) JacobiDN((alpha + eta1) TT (Xi - 2*(eta1^2 + alpha^2)) + 
       eK1, m1);
   Qminus(Xi_, 
     TT_) = ( (Gamma)sq - 1)/((Gamma)sq + 1) * (alpha + 
        eta1)/(alpha - eta1)* dnSlow*
     JacobiDN(2 eK1 (Abel - 1/4 + Omega(xi, TT)/(2 Pi)), m1);
   QminusPrime(Xi_, 
     TT_) = ( 2 (Gamma)sq)/((Gamma)sq + 1)^2 * (alpha + 
        eta1)*(kappa^2 + alpha *eta1)/R^2*dnSlow*
      JacobiDN(2 eK1 (Abel - 1/4 + Omega(Xi, TT)/(2 Pi)), m1)
     - 2 alpha *
      eta1/((alpha - eta1) R) *((Gamma)sq - 1)/((Gamma)sq + 1)*( 
       dnSlow*JacobiSN(2 eK1 (Abel - 1/4 + Omega(Xi, TT)/(2 Pi)), m1)*
         JacobiCN(2 eK1 (Abel - 1/4 + Omega(Xi, TT)/(2 Pi)), m1) + 
        snSlow*cnSlow*
         JacobiDN(2 eK1 (Abel - 1/4 + Omega(Xi, TT)/(2 Pi)), m1)
       );
   w22sq(Xi_, TT_) =  ((Gamma)sq + 1)^2/(4 (Gamma)sq) * 
     EllipticTheta(3, 0 , nome)^2*
     EllipticTheta(3, Pi (Abel + 1/4) + Omega(Xi, TT)/2 , 
        nome)^2/(EllipticTheta(3, Omega(Xi, TT)/2 , nome)^2*
        EllipticTheta(3, Pi (Abel + 1/4) , nome)^2);
   Y(Xi_, TT_) = (1 + Qminus(Xi, TT)^2)/(2 kappa);
   X(Xi_, TT_) = 
    1/(w22sq(Xi, TT)*(Sigma)*Exp(2*(CurlyPhi)(Xi, TT))) + 
     QminusPrime(Xi, TT);
   qsol(Xi_, 
     TT_) = (2 (1 - Qminus(Xi, TT)^2) X(Xi, TT) + 
       4 Qminus(Xi, TT)* Y(Xi, TT))/(X(Xi, TT)^2 + Y(Xi, TT)^2);
   
   qsol(xi, tval) + qBackground(xi, tval) 
   )
  )

SolSpeed(xi_?NumericQ, eta1_, eta2_, kappa_) :=
  Module({alpha},
   alpha = ComputeAlpha(x$val/t$val, eta1, eta2);
   2 (eta1^2 + alpha^2) + 
    4 kappa^2 EllipticK(eta1^2/alpha^2) / 
      EllipticPi(eta1^2/kappa^2, eta1^2/alpha^2) 
   );


SolPeakLeading(eta1_, eta2_, kappa_, x0_, sigma_, t$right_) :=
 
 Module({t$initial, x$initial, t$left, ODEeqns, xmax},
  t$initial = (-x0 + 3)/(4 kappa^2 - 4 eta1^2);
  x$initial = 
   NArgMax(sol(xmax, t$initial, eta1, eta2, kappa, x0, sigma), xmax);
  t$left = -x0/(4 kappa^2 - 4 eta1^2) + .1;
  ODEeqns =
   {
xsol'(tsol) == SolSpeed(xsol(tsol)/tsol, eta1, eta2, kappa),
x(t$initial) == x$initial
};
  NDSolve(ODEeqns, xsol, {tsol, t$left, t$right})
  )



SolPeakLeading(.25, 1, 2, -64, 1, 30)

NDSolve::ndode: The equations {True,0==SolSpeed(xsol(tsol)/tsol,0.25,1,2)} are not differential equations or initial conditions in the dependent variables {xsol}.
```