(*
Solution to Problem 1(e) of the Bigger Problem Set - KdV equation.
At the Mathematica prompt, type
<< kdv.math
*)

A[a_] := 4/3 a^2;

(* This is the analytic solution *)
da[th_] := A[a]/Cosh[a th]^2;

(* Check the analytic solution satisfies the KdV equation; chk should = 0 *)
chk = Simplify[ (- A[a]/2 + 3/2 da[th]) da'[th] + da'''[th]/6 ];

(* The differential equation together with its boundary conditions at th = 0 *)
deq[a_] := {(- A[a]/2 + 3/2 d[th]) d'[th] + d'''[th]/6 == 0, \
	d[0] == A[a], d'[0] == 0, d''[0] == - 2 A[a] a^2};

(* Solve the diffeq for the case a = 1 *)
soln = NDSolve[deq[1], d, {th,0,6}];

(* Plot the (log of the) numerical and analytic solutions *)
Plot[{Log[d[th] /. soln] , Log[da[th] /. a -> 1]}, {th,0,6}]