% This produces a black waterfall plot of the solution of 
% u_t + 6u^2u_x + u_xxx - (b(x,t)u)_x = 0 on [-pi,pi]
% and returns the the time vector to use for efft.m 
% The code is a simple modification of Trefethen's example p27.m
% The initial data is given by a double soliton with parameters
% a1, a2 (centers) and c1, c2 (velocities).
% The arguments of mkdvb.m are time t, parameters given in a vector form
% [a1,a2,c1,c2], a function bb = b(x,t), and the log_2 of the number of 
% x grid points (the default value is 8).
% It return the time vector T, the position vector X, and the matrix of
% solution values at times given by T.
% 
% Here is an example
% 
%  T = mkdvB(@(x,t) 100*cos(x-10^3*t).^2-50*sin(2*x+10^3*t), 0.05, [1,-1,4,-11],9)
%
% The time vector T can then be entered into effdyn.m code to compare the
% solution to the effective dynamics:
%
% [T,Y]=effdyn(@(x,t) 100*cos(x-10^3*t).^2-50*sin(2*x+10^3*t), T, [1,-1,4,-11])
%

function [T,X,U] = mkdvb2(fun,tmax,ZZ,N)

     if (nargin < 2 ) 

     tmax=0.065;

     end

     if (nargin < 3)

     ZZ = [-1,-2,-6,9];

     end

     if (nargin < 4)

     N = 8;

     end

N = 2^N;

% Set up grid and two-soliton initial data:

dt = .4/N^2; x = (2*pi/N)*(-N/2:N/2-1)';

drawnow, set(gcf,'renderer','zbuffer')

u=sqrt(2)*qu(x,ZZ);

v = fft(u); k = [0:N/2-1 0 -N/2+1:-1]'; ik3 = 1i*k.^3;

% Solve PDE and plot results:
nplt = floor((tmax/25)/dt); nmax = round(tmax/dt);
udata = u; tdata = 0; h = waitbar(0,'please wait...');
for n = 1:nmax
t = n*dt; g = -i*dt*k; 

bb=fun(x,t);

E = exp(dt*ik3/2); E2 = E.^2;
a = g.*fft(real( ifft(     v    ) ).^3 - ifft( v ).*bb);
b = g.*fft(real( ifft(E.*(v+a/2)) ).^3 - ifft(E.*(v+a/2)).*bb); 
c = g.*fft(real( ifft(E.*v + b/2) ).^3 - ifft(E.*v + b/2).*bb);     
d = g.*fft(real( ifft(E2.*v+E.*c) ).^3 -  ifft(E2.*v+E.*c).*bb);
v = E2.*v + (E2.*a + 2*E.*(b+c) + d)/6;

if mod(n,nplt) == 0 
u = real(ifft(v)); waitbar(n/nmax)
     udata = [udata u]; tdata = [tdata t];
    end
  end

udata = udata/sqrt(2); 

hold on 
     [K,P] = size(tdata);

z1=0;
z2=0;

for j=1:P
    plot3(x,tdata(j)+0*x, udata(:,j) );
z1=max(z1,max(udata(:,j)));
z2=min(z1,min(udata(:,j)));
end
T = tdata;
U = udata;
X = x;

view(0,60)

xlabel x, ylabel t, axis([-pi pi 0 tmax z2 z1]), grid off
close(h)

function qu = qu(x,ZZ)

     k1 = ZZ(3)*(x - ZZ(1));
     k2 = ZZ(4)*(x - ZZ(2));

     ga1 =  exp(-k1);
     ga2 = - exp(-k2);

     m11=(1+ga1.^2)./(2*ZZ(3));
     m22=(1+ga2.^2)./(2*ZZ(4));
     m12=(1+ga1.*ga2)./(ZZ(3)+ZZ(4));


     M=m11.*m22 - m12.^2;
     M1=ga1.*(m12 - m22) + ga2.*(m12-m11);
     qu = M1./M;

