# Unexpected stiff system error [closed]

Let there be the following functions

ya[t_, x_] := Cosh[t*x]
za[t_, x_] := Sinh[t*x]
wa[t_, x_] := x*Sinh[t*x]^2


One can verify that the following PDE's hold

D[ya[t, x], t] == wa[t, x]/za[t, x]

D[za[t, x], t] ==
D[wa[t, x], x]/
D[ya[t, x], x] + (x Cosh[t x] - (
Csch[t x] (2 t x Cosh[t x] Sinh[t x] + Sinh[t x]^2))/t)

D[wa[t, x], t] ==
D[wa[t, x], x]/
za[t, x] + (2 x^2 Cosh[t x] Sinh[t x] -
Csch[t x] (2 t x Cosh[t x] Sinh[t x] + Sinh[t x]^2))


i.e. that they are evaluated to True.

I would expect that the following NDSolve code would give ya, za and wa as an output:

mol[n_Integer, o_: "Pseudospectral"] := {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n,
"MinPoints" -> n, "DifferenceOrder" -> o}}

Monitor[AbsoluteTiming[

s = NDSolve[{

D[y[t, x], t] == w[t, x]/z[t, x],
D[z[t, x], t] ==
D[w[t, x], x]/
D[y[t, x], x] + (x Cosh[t x] - (
Csch[t x] (2 t x Cosh[t x] Sinh[t x] + Sinh[t x]^2))/t),
D[w[t, x], t] ==
D[w[t, x], x]/
z[t, x] + (2 x^2 Cosh[t x] Sinh[t x] -
Csch[t x] (2 t x Cosh[t x] Sinh[t x] + Sinh[t x]^2)),

w[t, 2] == wa[t, 2], y[t, 2] == ya[t, 2], z[t, 1] == za[t, 1],
w[1, x] == wa[1, x], y[1, x] == ya[1, x], z[1, x] == za[1, x]

}, {y, z, w}, {x, 1, 2}, {t, 1, 2}, Method -> mol[81],

EvaluationMonitor :> (currentTime =
Row[{" t = ", CForm[t]}])];], currentTime]


However an error occurs:

NDSolve::ndsz: At t == 1.071045212039208, step size is effectively zero; singularity or stiff system suspected.

I plotted ya, za and wa and they did not seem to be stiff at all, at least not at t == 1.071045212039208. So I think that there must be some kind of spelling mistake in my code.

I have not yet been able to find out what is going on. This is a real problem since it is this same error that I keep getting in other circumstances where an analytical solution is not available.

I would appreciate any help.

• The problem seems to be in the method you chose, not in the equations. If you remove the Method option from NDSolve and let it choose for itself, it finds a solution in a fraction of the time and with no errors. A plot of the interpolation function obtained for y is practically indistinguishable from ya, so the solution seems valid too. See: Plot3D[Evaluate[{y[t, x], ya[t, x]} /. First@s], {t, 1, 2}, {x, 1, 2}]. – MarcoB Jan 14 at 23:26
• I tried without Method but got NDSolve::ndnum: Encountered non-numerical value for a derivative at t == 1.. – jheidk51 Jan 15 at 11:13
• Is it possible that you had not defined ya, wa, and za when you ran this last try? – MarcoB Jan 15 at 13:41
• Yes it is. Sorry.. – jheidk51 Jan 15 at 13:53