function kdvFD_periodic
% Code written by D.N.Sibley (Loughborough University) for the SIG meeting
% of "Mathematical Challenges of Nonlinear Waves and Interfacial Dynamics",
% on 25th and 26th September 2017.
%
% Solves 1D KdV using central finite differences (with N points) subject to
% periodic BCs. tfinal is the end time point. Domain is x=-a to x=a.
N = 300; a = 20; tfinal = 1;
% So the stepsize and the x spatial domain are
dx = 2*a/(N-1); x = -a:dx:a;

% Set up the first and third order differentiation matrices using central
% differences, with periodic boundary conditions.
D1 = full(gallery('tridiag',N,-1,0,1));
D1(N,1)=1; D1(1,N)=-1;
D1 = D1/(2*dx);
% The 'tridiag' matrix doesn't have a natural extension to a quindiagonal
% form, so repeated use of diag implemented instead.
D3 =  diag( ones(N-2,1) , 2 )   - diag( ones(N-2,1) , -2 )...
    - diag( 2*ones(N-1,1) , 1 ) + diag( 2*ones(N-1,1) , -1 );
D3(1,N) = 2;  D3(1,N-1) = -1;
D3(2,N) = -1; D3(N-1,1) = 1;
D3(N,1) = -2; D3(N,2) = 1;
D3 = D3/(2*dx^3);

% Set the initial condition
u_initkdv = 6*sech(x).^2;

kdv_sol = ode15s(@kdv_eq,[0 tfinal],u_initkdv);

times = linspace(0,tfinal,100);
ukdv = deval(kdv_sol,times);

figure; waterfall(x,times,ukdv'); colormap([0 0 1]);
xlabel('$x$','FontSize',18,'interpreter','latex')
ylabel('$t$','FontSize',18,'interpreter','latex')
zlabel('$u(x,t)$','FontSize',18,'interpreter','latex')
set(gca,'fontsize',17)

pause; figure;

for i = 1:length(times)
    plot(x,ukdv(:,i),'b','LineWidth',2)
    hold on
    plot(x,12*(3 + 4*cosh(2*x-8*times(i))+cosh(4*x-64*times(i)))...
    ./(3*cosh(x-28*times(i))+cosh(3*x-36*times(i))).^2,'r:','LineWidth',2);
    hold off
    str=['KdV solitons ','$t = \,$',num2str(times(i),'%1.2f')];
    title(str,'interpreter','latex')
    xlabel('$x$','FontSize',18,'interpreter','latex')
    ylabel('$u(x)$','FontSize',18,'interpreter','latex')
    set(gca,'fontsize',17)
    axis([-a a -5 10])
    legend('Numeric','Analytic','location','northwest');
    pause(0.01)
end

% The KdV eqn
function udot = kdv_eq(~,u)
    udot = - 6*u.*(D1*u) - (D3*u);
end

end

% Matlab commands to read using "help ***" (if interested):
%
% full; gallery; diag; ones; ode15s; linspace; deval; figure; waterfall;
% colormap; xlabel; set; pause; length; plot; hold; num2str; title; axis;
% legend;
%