function [ x, sol, sig, ta ] = cburgermol( n )
%CBURGERMOL Solve complex Burger's eqn
%   u_t + u u_x = 0
%   -1 <= x <= 1, u(t,-1) = u(t,1)
%   u(0,x) = exp( i*pi*x )
%   JAC Weideman's solution u(t,x) = W( i*pi*t*exp(i*pi*x) )/(i*pi*t)
%   t > 0, W = LambertW.  Singular at x=1/2, t=1/pi/e.
%
%   u = sigma + i*tau gives
%   
%  [ sigma_t ] = - [ sigma  -tau  ][ sigma_x ]
%  [  tau_t  ]     [  tau   sigma ][  tau_x  ]
x = cos( pi*(0:n)/n );
D = gallery('chebspec', n+1 );
inif = [ cos(pi*x), sin(pi*x) ];
% Spectral space derivatives
% Enforce periodic boundary conditions by using the mean
% of the left and right end-values to advance the left and 
% right end-values.
    function yp = cb(t,y)
        yp = zeros(size(y));
        sigma = y(1:n+1,:);
        tau   = y(n+2:end,:);
        sp = (sigma(1,:)+sigma(n+1,:))/2;
        sigma(1,:)=sp;
        sigma(n+1,:)=sp;
        tp = (tau(1,:)+tau(n+1,:))/2;
        tau(1,:)=tp;
        tau(n+1,:)=tp;
        sx = D*sigma;
        tx = D*tau;
        yp(1:n+1,:) = -(sigma.*sx - tau.*tx);
        yp(n+2:end,:) = -(tau.*sx + sigma.*tx);
    end

opts = odeset('reltol',1.0e-6);
ts = 1/pi/exp(1);
sol = ode15s(@cb,[0,ts], inif, opts );

% Postprocess solution for visualization
sig = sol.y(1:n+1,:);
ta = sol.y(n+2:end,:);
figure
mesh(x,sol.x,sig');
xlabel('x','fontsize',16);
ylabel('t','fontsize',16);
zlabel('real','fontsize',16);
set(gca,'fontsize',16);
axis([-1,1,0,ts,-1,1]);
view(170,40);
figure
mesh(x,sol.x,ta');
xlabel('x','fontsize',16);
ylabel('t','fontsize',16);
zlabel('imag','fontsize',16);
set(gca,'fontsize',16);
axis([-1,1,0,ts,-1,3]);
end

