3 Added an update on Mathematica 10 and an animated GIF
source | link

Update 2
After upgrading to version 10, I discovered that Mathematica was trying to be too careful. Well, the support staff at Wolfram told me, but anyways the code needed to be altered to prevent a memory overflow. The simple solution was to do the following:

eqall = N[{eq0, eq1, eq2, eq3, eq4a, eq4b}];

Apparently, using what Mathematica defines as exact numbers caused a problem. Mathematica would finish the computation, but the memory would continue to disappear. After the addition of N[], the issue went away. Interestingly, at the point where the wave started to turn over, the code now fails. It had always ran into stiff solutions, as it should. But before it did not hit this limit until t ~ 51, whereas now it runs into problems at t > 47.95. In any case, I made a GIF that shows the steepening evolve, which is what I was going for in the first place.

enter image description here

Update 2
After upgrading to version 10, I discovered that Mathematica was trying to be too careful. Well, the support staff at Wolfram told me, but anyways the code needed to be altered to prevent a memory overflow. The simple solution was to do the following:

eqall = N[{eq0, eq1, eq2, eq3, eq4a, eq4b}];

Apparently, using what Mathematica defines as exact numbers caused a problem. Mathematica would finish the computation, but the memory would continue to disappear. After the addition of N[], the issue went away. Interestingly, at the point where the wave started to turn over, the code now fails. It had always ran into stiff solutions, as it should. But before it did not hit this limit until t ~ 51, whereas now it runs into problems at t > 47.95. In any case, I made a GIF that shows the steepening evolve, which is what I was going for in the first place.

enter image description here

2 update about version compatibility issues with, possibly, new java versions
source | link

Update
After I updated java to Version 7 Update 71, the above code stopped working (i.e., it would quit the kernel and undefine everthing as if I had just opened the notebook) for Mathematica v8.0.1 on both my laptop (running OS X 10.9.5) and my desktop (running OS X 10.8.5). I updated my laptop to Mathematica 10 over the weekend and tried again. Then the Mathematica code worked again. My only guess as to the source of the problem was a change in the garbage collection of java, but JLink said that I was still using an earlier version of java. I wrote a formal report to Wolfram last Monday, but never heard back so I "bit the bullet" and upgraded hoping that would solve the problem. It did work, but I would have appreciated knowing this 5 days earlier so I would not have wasted a week I did not have trying to fix the issue.

I am not sure if the java update caused the problems, but I could not see why the bash update would have any effect on Mathematica. Then again, I am not an expert, so take my comments with a grain of salt.

Update
After I updated java to Version 7 Update 71, the above code stopped working (i.e., it would quit the kernel and undefine everthing as if I had just opened the notebook) for Mathematica v8.0.1 on both my laptop (running OS X 10.9.5) and my desktop (running OS X 10.8.5). I updated my laptop to Mathematica 10 over the weekend and tried again. Then the Mathematica code worked again. My only guess as to the source of the problem was a change in the garbage collection of java, but JLink said that I was still using an earlier version of java. I wrote a formal report to Wolfram last Monday, but never heard back so I "bit the bullet" and upgraded hoping that would solve the problem. It did work, but I would have appreciated knowing this 5 days earlier so I would not have wasted a week I did not have trying to fix the issue.

I am not sure if the java update caused the problems, but I could not see why the bash update would have any effect on Mathematica. Then again, I am not an expert, so take my comments with a grain of salt.

1
source | link

This question relates to nonlinear wave steepening and I see now that I did not list everything that was necessary to address the problem. I am also quite certain that I phrased this question poorly on account of my novice-level experience with NDSolve in Mathematica.

The question spawned from my attempt to reproduce results from this paper. On the second page, they go into some details on how they analyzed the problem by simplifying their original quasi-general magnetized fluid equations. On the right-hand side of this page, the authors seemed to imply that they produced the results in their Figure 1 using the approach I outlined in my OP. Unfortunately, in my frustration I failed to catch an important detail on the last page of their paper.

The authors analyzed two equations numerically, shown here: $$ \begin{align} \left( \partial{\scriptstyle_{\tau}} + w \ \partial{\scriptstyle_{\chi}} \right) w + \partial{\scriptstyle_{\chi}} \left( h + \beta \ln h \right) & = 0 \\ \partial{\scriptstyle_{\tau}} \ h + h \ \partial{\scriptstyle_{\chi}} \ w + w \ \partial{\scriptstyle_{\chi}} \ h & = 0 \end{align} $$ where $h$ represents a normalized magnetic field component and $w$ a normalized velocity component. They assumed the following boundary conditions: $$ \begin{align} \lim_{\lvert x \rvert \rightarrow \infty} w\left( \chi, \tau \right) & = 0 \\ \lim_{\lvert x \rvert \rightarrow \infty} h\left( w\left( \chi, \tau \right) \right) & = 1 \\ \end{align} $$ and assumed the following initial conditions: $$ \lim_{\tau \rightarrow 0} h\left( w\left( \chi, \tau \right) \right) = 1 + \frac{ 1 }{ 2 } \cosh^{-1}{ \left( \frac{ \chi }{ 20 } \right) } $$ I finally realized, after wasting a considerable amount of time, that the pieces I was missing were: the first two PDEs in this answer and the initial condition for $w(\chi,0)$. The last missing piece is given by: $$ w = 2 \left( \sqrt{ h + \beta } - \sqrt{ 1 + \beta } \right) + \sqrt{ \beta } \left( \ln \lvert \frac{ \sqrt{ h + \beta } - \sqrt{ \beta } }{ \sqrt{ h + \beta } + \sqrt{ \beta } } \rvert - \ln \lvert \frac{ \sqrt{ 1 + \beta } - \sqrt{ \beta } }{ \sqrt{ 1 + \beta } + \sqrt{ \beta } } \rvert \right) $$ where we substitute our initial condition for $h(\chi,0)$ shown above into this equation. Then, we need only let Mathematica do the work. The code below works, though it is computationally expensive and made my computer work enough to get the fans going.

tmax = 100.;  (* Max. time step *)
xmax = 500.;  (* Max. spatial point *)
dtmaxf = 0.0001;  (* Max. fractional step *)
beta = 2.;  (* value for plasma beta *)
xfac = 20; (* normalization factor for Sech *)
xbounds = {-xmax/2, xmax};
tbounds = {-1, tmax};
(* Define initial condition functions *)
hatt0[x_, xo_] := 1 + 1/2 Sech[x/xo];
wofhb[hh_, bb_] := Sqrt[bb] + Sqrt[1 + bb] (-2 Sqrt[1 + bb] + 2 Sqrt[hh + bb] -
                   Sqrt[bb] Log[(-Sqrt[bb] + Sqrt[1 + bb])/()] + Sqrt[bb]
                   Log[(-Sqrt[bb] + Sqrt[hh + bb])/(Sqrt[bb] + Sqrt[hh + bb])]);
(* Define partial differential equations *)
eq0 = D[w[x, t], t] + w[x, t] D[w[x, t], x] + D[(h[x, t] + be Log[h[x, t]]), x] == 0;
eq1 = D[h[x, t], t] + w[x, t] D[h[x, t], x] + h[x, t] D[w[x, t], x] == 0;
eq2 = h[x, 0] == hatt0[x, xfac];
eq3 = w[x, 0] == (wofhb[h, be]) /. {h -> (hatt0[x, xfac])};
eq4a = h[xmax, t] == 1;
eq4b = w[xmax, t] == 0;
eqall = {eq0, eq1, eq2, eq3, eq4a, eq4b};

(* Find solutions to partial differential equations as function of [x,t] *)
soln = Block[{c = 0}, NDSolve[(eqall /. {be -> beta}),
             {w[x, t], h[x, t]},
             {x, xbounds[[1]], xbounds[[2]]},
             {t, tbounds[[1]], tbounds[[2]]},
             MaxSteps -> 10^6, SolveDelayed -> True, 
             Method -> {"MethodOfLines", "TemporalVariable" -> t}, 
             MaxStepFraction -> dtmaxf, SolveDelayed -> True, 
             StepMonitor :> If[Mod[c, 1000] == 0, Print["Step # = ", c]; c++, c++]]]

Then we can plot these results using the following:

(* Define time steps and colors for plots *)
tsolns = {0, 10, 20, 30, 40, 50};
cols = {Black, Purple, Blue, Green, Orange, Red};
(* Define plot outputs *)
plots = Table[Plot[((Evaluate[(h[x, t] /. Flatten[soln][[2]])]) /. {t -> tsolns[[i]]}), 
              {x, -100, 180}, PlotRange -> All, 
              PerformanceGoal -> "Quality", 
              PlotStyle -> {AbsoluteThickness[1], cols[[i]]}], 
              {i, 1, Length[tsolns], 1}];
(* Show plots *)
Show[plots, ImageSize -> {600, 600}, ImageMargins -> 0, AspectRatio -> Full]

The results are shown in the image below. Note that the last (red) time step shows a "fuzzy" part near x ~ 150. This is an error (perhaps not the right word) due to the steepening reaching a gradient catastrophe. At that point, the results become multivalued and I am not clever enough with Mathematica to know how to deal with this. I assume it has to do with the Stiffness settings in NDSolve, but again, I was not sure. I was just happy to get this working as well as it did. If any of the Mathematica whisperers on this site have a more efficient suggestion and a better way to handle the gradient catastrophe, I would greatly appreciate hearing about it.

Plot Description
The plot shows six snapshots of $h(\chi,\tau)$ at $\tau$ = 0 (black), 10 (purple), 20 (blue), 30 (green), 40 (orange), and 50 (red) versus $\chi$.

Snapshots of steepening wave