function DoThemAll ()
% This is the Method of Lines for the EW equation.
% Version without explicit boundary conditions. 
% We use the starting function as a trial function.
% Written by R. Schaback, Feb. 2010
clear all; % delete garbage from previous runs
close all; % close old figures, if they exist
% Some globals, see EWRSnbd.m below
global c;   % scale of kernel
global U    % U of paper
global Up   % U prime of paper
global Bdagger % the pseudoinverse, see the paper
global u0 % starting function
global u0p % first derivative thereof
% Some other globals concerning RBFs,
% see frbf.m below
global RBFscale
global RBFtype
global RBFpar
RBFtype='g'; % Gaussian is 'g', others see frbf.m below
RBFpar=-2;
% standard parameters of equal width PDE
mu=1/6;
k=sqrt(1/(4*mu));
userinput=input('Enter single problem number between 1 and 7, or 0 for all:>');
movieinput=input('Enter 1 if you want to see movies, or 0 for none:>');
timebreaks=51; % breaks for plotting when ODE solving, should be at least 2
% Note that, if movies are requested, this is the number of figures
% generated
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now set up the problem 
%
% Define a,b, and endtime. Starttime is always 0.
% ptr must be ntr=n tRial points in space
% pts will be nts=N>=n  teSt points in space, also used for plotting
% We need u0 and its first derivative u0p
% at the test points pts. 
% maxy should ne roughly the max of what is to be plotted.
% miny should ne roughly the max of what is to be plotted.
% plotpos controls where the legend appears.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for icase=1:6
    if userinput~=0 
        icase=userinput;
    end
    switch (icase)
        case 4
            % undular bore
            a=-20;
            b=50;
            endtime=800;
            ntr=200;
            ptr=(a:(b-a)/(ntr-1):b)'; % spatial trial points, n=ntr in paper
            nts=500;
            disp(sprintf('Undular bore on [%d, %d] until T=%d on %dx%d points',...
                a, b, endtime, ntr, nts));  
            pts=(a:(b-a)/(nts-1):b)'; % spatial interpolation points, N=nts in paper
            xc=0;
            d=2;
            maxy=0.2;
            miny=-0.1;
            u00=0.1;
            u0=(u00/2)*(1-tanh((pts-xc)/d)); % starting function and derivatives
            u0p =-1/2*u00*(1-tanh((pts-xc)/d).^2)/d;
            plotpos=3;
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % Gaussian (Maxwellian) starting function
        case 5
            a=0;
            b=50;
            endtime=500;
            xc=20
            d=5;
            ntr=200;
            ptr=(a:(b-a)/(ntr-1):b)'; % spatial trial points, n=ntr in paper
            nts=500;
            disp(sprintf('Maxwellian case  on [%d, %d] until T=%d on %dx%d points',...
                a, b, endtime, ntr, nts));  
            pts=(a:(b-a)/(nts-1):b)'; % spatial interpolation points, N=nts in paper
            maxy=0.2;
            miny=-0.2;
            u00=0.1;
            u0 =1/2*u00*exp(-(pts-xc).^2/d^2);
            u0p=-u00*(pts-xc)/d^2.*exp(-(pts-xc).^2/d^2);
            plotpos=1;
        case 1
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % Single soliton starting function, zero boundary values, Raslan paper
            a=0;
            b=30;
            endtime=80;
            xc=10;
            mu=1;
            p=0.5/sqrt(mu);
            cc=0.1;
            maxy=0.4;
            miny=-0.05;
            ntr=200;
            ptr=(a:(b-a)/(ntr-1):b)'; % spatial trial points, n=ntr in paper
            nts=500;
            disp(sprintf('Single soliton on [%d, %d] until T=%d on %dx%d points',...
                a, b, endtime, ntr, nts));  
            pts=(a:(b-a)/(nts-1):b)'; % spatial interpolation points, N=nts in paper
            u0 = 3* cc*sech(p*(pts-xc)).^2;
            u0p=-6* cc*sech(p*(pts-xc)).^2.*tanh(p*(pts-xc))*p;
            plotpos=1;
            
        case 7
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % Double soliton starting function, 
            % positive and negative, zero boundary values, Raslan paper
            a=-10;
            b=70;
            endtime=30;
            xc=23;
            mu=1;
            p=0.5/sqrt(mu);
            cc=1.2;
            maxy=6;
            miny=-6;
            ntr=200;
            ptr=(a:(b-a)/(ntr-1):b)'; % spatial trial points, n=ntr in paper
            nts=500;
            disp(sprintf('Double clashing solitons on [%d, %d] until T=%d on %dx%d points',...
                a, b, endtime, ntr, nts));  
            
            pts=(a:(b-a)/(nts-1):b)'; % spatial interpolation points, N=nts in paper
            u01 = 3* cc*sech(p*(pts-xc)).^2;
            u01p=-6* cc*sech(p*(pts-xc)).^2.*tanh(p*(pts-xc))*p;
            cc=-1.4;
            xc=38;
            u02 = 3* cc*sech(p*(pts-xc)).^2;
            u02p=-6* cc*sech(p*(pts-xc)).^2.*tanh(p*(pts-xc))*p;           
            u0=u01+u02;
            u0p=u01p+u02p;
            plotpos=1;
            
        case 2
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % Double soliton starting function, 
            % positive, zero boundary values, Raslan paper
            a=-10;
            b=70;
            endtime=30;
            xc=10;
            mu=1;
            p=0.5/sqrt(mu);
            cc=1.5;
            maxy=6;
            miny=-0.05;
            ntr=200;
            ptr=(a:(b-a)/(ntr-1):b)'; % spatial trial points, n=ntr in paper
            nts=500;
            disp(sprintf('Double solitons on [%d, %d] until T=%d on %dx%d points',...
                a, b, endtime, ntr, nts));  
            
            pts=(a:(b-a)/(nts-1):b)'; % spatial interpolation points, N=nts in paper
            u01 = 3* cc*sech(p*(pts-xc)).^2;
            u01p=-6* cc*sech(p*(pts-xc)).^2.*tanh(p*(pts-xc))*p;
            u01pp=12*cc*sech(p*(pts-xc)).^2.*tanh(p*(pts-xc)).^2*p^2-6*cc*sech(p*(pts-xc)).^2.*(1-tanh(p*(pts-xc)).^2)*p^2;
            cc=0.75;
            xc=25;
            u02 = 3* cc*sech(p*(pts-xc)).^2;
            u02p=-6* cc*sech(p*(pts-xc)).^2.*tanh(p*(pts-xc))*p;
            u0=u01+u02;
            u0p=u01p+u02p;
            plotpos=1;
            
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        case 3
            % Triple soliton starting function, 
            % positive, zero boundary values, Raslan paper
            a=-10;
            b=100;
            endtime=15;
            xc=10;
            mu=1;
            p=0.5/sqrt(mu);
            cc=4.5;
            maxy=15;
            miny=-0.2;
            ntr=200;
            ptr=(a:(b-a)/(ntr-1):b)'; % spatial trial points, n=ntr in paper
            nts=500;
            disp(sprintf('Triple solitons on [%d, %d] until T=%d on %dx%d points',...
                a, b, endtime, ntr, nts));  
            
            pts=(a:(b-a)/(nts-1):b)'; % spatial interpolation points, N=nts in paper
            u01 = 3* cc*sech(p*(pts-xc)).^2;
            u01p=-6* cc*sech(p*(pts-xc)).^2.*tanh(p*(pts-xc))*p;
            cc=1.5;
            xc=25;
            u02 = 3* cc*sech(p*(pts-xc)).^2;
            u02p=-6* cc*sech(p*(pts-xc)).^2.*tanh(p*(pts-xc))*p;
            cc=0.5;
            xc=35;
            u03 = 3* cc*sech(p*(pts-xc)).^2;
            u03p=-6* cc*sech(p*(pts-xc)).^2.*tanh(p*(pts-xc))*p;
            u0=u01+u02+u03;
            u0p=u01p+u02p+u03p;
            plotpos=1;
            
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        case 6
            % Double soliton starting function, 
            % positive and negative, zero boundary values, RS, full kill
            a=-40;
            b=40;
            endtime=100;
            xc=20;
            mu=1;
            p=0.5/sqrt(mu);
            cc=-1.2;
            maxy=6;
            miny=-6;
            ntr=1001;
            ptr=(a:(b-a)/(ntr-1):b)'; % spatial trial points, n=ntr in paper
            nts=1001;
            disp(sprintf('Double clashing solitons on [%d, %d] until T=%d on %dx%d points',...
                a, b, endtime, ntr, nts));  
            
            pts=(a:(b-a)/(nts-1):b)'; % spatial interpolation points, N=nts in paper
            u01 = 3* cc*sech(p*(pts-xc)).^2;
            u01p=-6* cc*sech(p*(pts-xc)).^2.*tanh(p*(pts-xc))*p;
            cc=1.2;
            xc=-20;
            u02 = 3* cc*sech(p*(pts-xc)).^2;
            u02p=-6* cc*sech(p*(pts-xc)).^2.*tanh(p*(pts-xc))*p;
            u0=u01+u02;
            u0p=u01p+u02p;
            plotpos=1;
            
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        otherwise
            error('Illegal problem number')
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % End of problem setup
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % We first want to find a good scale for RBF trial functions 
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    c=3*(b-a)/ntr; % first guess for scale, good for Gaussian
    RBFscale=c;
    if RBFtype=='w' % for Wendland functions, use larger startup scales
        RBFscale=RBFscale*10;
    end
    % The scale is critical. We loop for it,
    % starting at large scales until condition is OK
    disp('Trying first to find a good scale for the kernel')
    while 1==1
        A=kermat(ptr(:,1),ptr(:,1)); % standard kernel matrix
        ce=condest(A); % condition estomate
        if ce<1.e7     % this condition level may be sensitive 
            break;
        else
            RBFscale=RBFscale/1.1;
        end
    end
    c=RBFscale;
    final_scale=c
    % This ends the scale estimation.
    disp('Preparing the spatial matrices')
    [U, Up, Upp]=U31Dkermat(pts(:,1),ptr(:,1));
    B=mu*Upp-U; % B like in the paper
    % We now get the pseudo-inverse
    Bdagger=pinv(B);
    % Now we prepare to solve the ODE system.
    startcoeff=zeros(ntr,1);
    disp('The spatial problem is set up now.')
    disp('Hit ENTER to start ODE solving')
    pause
    disp('Starting the ODE solver')
    % We now invoke an ODE solver,
    % a stiff one or a fast one....
    % [t,y] = ode15s(@EWRSnbd,0:endtime/(timebreaks-1):endtime,startcoeff);
    [t,y] = ode113(@EWRSnbd,0:endtime/(timebreaks-1):endtime,startcoeff);
    disp('ODE system is solved. Now for graphical output.')
    if movieinput==1 % now do a movie
        for i=1:length(t)
            figure
            plot(pts(:,1),u0+U*y(i,:)',pts(:,1),u0);
            axis([a b miny maxy])    
            F(i)=getframe;
            hold on
        end 
        figure
        axis([a b 0 0.2])    
        movie(F)
    end
    % now just the final black-and-white plot for the paper
    i=length(t);
    figure
    plot(pts(:,1),u0+U*y(i,:)','k-',pts(:,1),u0,'k:');
    axis([a b miny maxy])    
    legend(sprintf('Solution at time %d\n',endtime),'Starting function',plotpos)
    xlabel('x')
%     % now for the invariants, up to a fixed factor
%     for i=1:length(t)
%         uval = (u0+U*y(i,:)');
%         uvalp= (u0p+Up*y(i,:)');
%         invar(i,1)=sum(uval);   
%         invar(i,3)=sum(uval.*uval.*uval);   
%         invar(i,2)=sum(uval.*uval+mu*uvalp.*uvalp);
%         invar(i,4)=uval(1);
%         invar(i,5)=uvalp(1);
%     end
%     % figure
%     % plot(invar(:,1))
%     % title('Invariant 1')
%     % figure
%     % plot(invar(:,2))
%     % title('Invariant 2')
%     % figure
%     % plot(invar(:,3))
%     % title('Invariant 3')
%     inv1=max(abs(invar(:,1)-mean(invar(:,1))))
%     inv2=max(abs(invar(:,2)-mean(invar(:,2))))
%     inv3=max(abs(invar(:,3)-mean(invar(:,3))))
    disp('Hit ENTER to close all figure windows')
    pause
    close all;
    if userinput~=0 
        break;
    end
    
end
return

function dydt = EWRSnbd(t,y)
% This is the rhs of the ODE system for the EW equation
% For the globals, see the main program
global U
global Up
global Bdagger
global u0
global u0p
dydt=Bdagger*((u0+U*y).*(u0p+Up*y));
return

function [U, Up, Upp]=U31Dkermat(X, Y)
% creates 3 kernel matrices
% for two 1D point sets X and Y.
% U is the standard kernel matrix, Up the x derivative,
% Upp the double x derivative.
% See frbf.m for further details.
global RBFscale
[nx dx]=size(X);
[ny dy]=size(Y);
if dx~=dy
    error('Unequal space dimension for U31Dkermat arguments');
end
U=frbf(distsqh(X,Y)/RBFscale^2,0);
f1mat=frbf(distsqh(X,Y)/RBFscale^2,1)/RBFscale^2;
f2mat=frbf(distsqh(X,Y)/RBFscale^2,2)/RBFscale^4;
dmat=(repmat(X(:,1),1,ny)-repmat(Y(:,1)',nx,1));
Up=f1mat.*dmat;
Upp=f1mat+f2mat.*dmat.*dmat;
return

function y =  frbf (s, k)   
% k-th derivative of standard rbf kernel in f form, i.e.
% as a function of s=r^2/2.
% To understand what is happening here, consult
% Programming Hints for Kernel-Based Methods   on RS's website
% NOTE: in the Equal Width Equation paper, kernels should
% ALWAYS decay towards infinity.
global RBFtype;   % 'g'=gaussian, 'mq' = multiquadric
global RBFpar;
switch lower(RBFtype)
    case ('g')   % 'g'=gaussian, 
        if mod(k,2)==0
            y=exp(-s);
        else
            y=-exp(-s);
        end
    case ('mq') % 'mq' = multiquadric, inverse or not...
        fac=1;
        ord=k;
        par=RBFpar;
        while ord>0
            ord=ord-1;
            fac=fac*par;
            par=par-2;
        end
        y=fac*(1+2*s).^(par/2);    
    case ('p')  % powers 
        fac=1;
        ord=k;
        par=RBFpar;
        while ord>0
            ord=ord-1;
            fac=fac*par;
            par=par-2;
        end
        y=fac*(2*s+eps).^(par/2);  
    case ('tp')  % thin-plate 
        fac=1;
        ord=k;
        par=RBFpar;
        su=0;
        while ord>0
            ord=ord-1;
            if ord==k-1
                su=1;
            else
                su=su*par+fac;
            end           
            fac=fac*par;
            par=par-2;
        end
        y=(2*s+eps).^(par/2);
        y=fac*y.*log(2*s+eps)/2 +su*y;    
    case ('ms')  % Matern/Sobolev
        y=(-1)^k*besselk(RBFpar-k,sqrt(2*s+eps)).*(sqrt(2*s+eps)).^(RBFpar-k);
    case ('w') % Wendland functions.
        % we use only those which are pos. def. in dimension at most 3.
        [coeff, expon]=wendcoeff(3+2*k, RBFpar-k);
        r=sqrt(2*s);
        ind=find(r<=1);
        u=zeros(size(ind));
        sp=ones(size(ind));
        sloc=r(ind);
        for i=1:length(coeff)
            u=u+coeff(i)*sp;
            sp=sp.*sloc;   
        end
        u=u.*(1-sloc).^expon;
        y=zeros(size(s));
        y(ind)=(-1)^k*u;
    otherwise
        error('RBF type not implemented')
end
return

function dst=distsqh(p, q)
% calculates a np*nq matrix of halved squares of point distances 
% with two point sets p and q of np and nq points each.
[np pdim]=size(p);
[nq qdim]=size(q);
if pdim~=qdim 
    error('point sets of unequal dimension')
end
if pdim==1
    dst=(repmat(p,1,nq)-repmat(q',np,1)).^2/2;
    return
end
dst=p*q';
cp=sum((p.*p)')/2; % squared norms of p points, as row, halved
cq=sum((q.*q)')/2; % squared norms of q points, as row, halved
dst=repmat(cp',1,nq)+repmat(cq,np,1)-dst; 
return

function mat=kermat(X, Y)
% creates kernel matrix
% for two point sets X and Y  
global RBFscale
[nx dx]=size(X);
[ny dy]=size(Y);
if dx~=dy
    error('Unequal space dimension for kermat arguments');
end
mat=frbf(distsqh(X,Y)/RBFscale^2,0);
return


function [coeff, expon]=wendcoeff(d,k)
% calculates coefficients and exponent for polynomial 
% part p_{d,k}(r) of Wendland functions 
% phi_{d,k}(r)=p_{d,k}(r)*(1-r)^expon
expon=floor(d/2)+k+1;
coeff=zeros(k+1,1);
coeff(1,1)=1;
for n=0:k-1
    coeff(n+2,1)=coeff(n+1,1)/(n+expon+2);
    for j=n+1:-1:2
        coeff(j,1)=(j*coeff(j+1,1)+coeff(j-1,1))/(expon+j); 
    end
    expon=expon+1;
    coeff(1,1)=coeff(2,1)/expon;
end


