# 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(tinitial) == xinitial }; NDSolve(ODEeqns, xsol, {tsol, tleft, 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}.
$$```$$
``````