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}.
```
```