function [ y ] = shocktube5tmol( n )
%SHOCKTUBETMOL Solve shock tube gas equation by TMOL
%   Use bvp4c and Backward Euler
%   u = [ rho; rho_x; m; m_x; E; E_x ]
%   u_t = -(u^2/2)_x + nu* v_x
close all

gam = 1.4;
% Auxiliary routines for bvp solver
    function y = inivals( x )
        k = find(x<=0);
        m = find(x>0);
        y = [ones(size(x(k))),0.125*ones(size(x(m)))
             zeros(size(x))
             zeros(size(x))
             zeros(size(x))
             2.5*ones(size(x(k))),0.25*ones(size(x(m)))
             zeros(size(x))];
    end
    function res = bc( ya, yb )
        res = [ya(1)-1
               ya(3)
               ya(5)-2.5
               yb(1)-0.125
               yb(3)
               yb(5)-0.25];
    end
inimesh = [-5:-1,-0.5,-0.25,-0.125,-0.02,-0.01,-0.001,-0.0001,0.0001,0.001,0.01,0.02, 0.125,0.25,0.5,1:5];

%
% Finite difference Jacobians were too susceptible to singularity.
% The analytical Jacobian below was computed in Maple and then 
% converted to Matlab via Maple's CodeGeneration[Matlab] facility.
%
% Nostalgically, this Jacobian was one of the very first uses of Maple,
% ever, and was significant for the 1982 thesis. RMC had forgotten that,
% and was pleased to see that Maple is still helpful for the same
% reasons today.  (And a little more so, too).
%
  function J = jacoshock( x, y )
        rho   = y(1,:);
        rhox  = y(2,:);
        m     = y(3,:);
        mx    = y(4,:);
        E     = y(5,:);
        Ex    = y(6,:);
        J = [0 1 0 0 0 0
            0.1e1 / h / nu 0 0 1 / nu 0 0
            0 0 0 1 0 0
            2*m*(-2*mx*rho+mx*rho*gam+2*m*rhox-m*gam*rhox)/rho^3/nu ...
            m^2*(-2+gam)/rho^2/nu ...
            -(-4*mx*rho*h+2*gam*mx*rho*h+4*m*rhox*h-2*m*gam*rhox*h-rho^2)/rho^2/h/nu ...
            -2*m*(-2+gam)/rho/nu 0 (gam-1)/nu
            0 0 0 0 0 1
            -(mx * gam * E * rho ^ 2 - 6 * mx * rho * gam * m ^ 2 + 6 * mx * rho * m ^ 2 ...
            + m * gam * Ex * rho ^ 2 - 2 * m * gam * E * rhox * rho + 6 * m ^ 3 * rhox * gam - 6 * m ^ 3 * rhox) / rho ^ 4 / nu ...
            -m * (gam * E * rho - 2 * m ^ 2 * gam + 2 * m ^ 2) / rho ^ 3 / nu ...
            (-6 * mx * m * rho * gam + 6 * mx * m * rho + gam * Ex * rho ^ 2 - rho * gam * E * rhox + 6 * rhox * gam * m ^ 2 - 6 * rhox * m ^ 2) / rho ^ 3 / nu ...
            (gam * E * rho - 3 * m ^ 2 * gam + 3 * m ^ 2) / rho ^ 2 / nu ...
            (gam * mx * rho * h - m * gam * rhox * h + rho ^ 2) / rho ^ 2 / h / nu m / rho * gam / nu];

    end



h = 5.0e-3; % As in Corless & Simpson (1983)

y = cell(n);
% Initial values duplicated in this function;
% if the trapezoidal rule is used, we need 2nd x derivative too.
    function yp = shocktube1( x, y )
        rhsx = inivals(x);
        yp = zeros(size(y));
        rho   = y(1,:);
        rho_x = y(2,:);
        m     = y(3,:);
        m_x   = y(4,:);
        E     = y(5,:);
        E_x   = y(6,:);
        yp(1,:) = rho_x ;
        yp(3,:) = m_x;
        yp(5,:) = E_x;
        p = (gam-1)*(E-(m.^2 ./rho )/2);
        p_x = (gam-1)*(E_x - m.*m_x./rho + m.^2 .*rho_x./rho.^2/2);
        yp(2,:) = (m_x)/nu + ( rho - rhsx(1,:) )/h/nu;
        yp(4,:) = (2*m.*m_x./rho-m.^2.*rho_x./rho.^2+p_x)/nu + ( m - rhsx(3,:) )/h/nu;
        yp(6,:) = (m_x.*(E+p)./rho + m.*(E_x+p_x)./rho - m.*(E+p).*rho_x./rho.^2)/nu + ( E - rhsx(5,:) )/h/nu;
    end
% Take the first step from the initial condition
opts = bvpset('FJacobian',@jacoshock,'reltol',1.0e-3);
nu = 8.0e-2;
solinit = bvpinit( inimesh, @inivals );
y{1} = bvp4c( @shocktube1, @bc, solinit, opts );
% Gradually reduce nu, because initial step with shock is too much
% for this method.
for i=1:4,
    nu = nu/2;
    h = h/4;
    solinit = bvpinit( y{1}, [-5,5] );
    y{1} = bvp5c( @shocktube1, @bc, solinit, opts );
end

% Subsequent steps use the previous solution
% as initial conditions, and also as the right-hand side in the BVP.
% This is passed to the following function by nested scoping.
    function yp = shocktube( x, y )
        yo = deval( rhssol, x );
        rhsx = yo;
        yp = zeros(size(y));
        rho   = y(1,:);
        rho_x = y(2,:);
        m     = y(3,:);
        m_x   = y(4,:);
        E     = y(5,:);
        E_x   = y(6,:);
        yp(1,:) = rho_x ;
        yp(3,:) = m_x;
        yp(5,:) = E_x;
        p = (gam-1)*(E-(m.^2 ./rho )/2);
        p_x = (gam-1)*(E_x - m.*m_x./rho + m.^2 .*rho_x./rho.^2/2);
        yp(2,:) = (m_x)/nu + ( rho - rhsx(1,:) )/h/nu;
        yp(4,:) = (2*m.*m_x./rho-m.^2.*rho_x./rho.^2+p_x)/nu + ( m - rhsx(3,:) )/h/nu;
        yp(6,:) = (m_x.*(E+p)./rho + m.*(E_x+p_x)./rho - m.*(E+p).*rho_x./rho.^2)/nu + ( E - rhsx(5,:) )/h/nu;
    end

% Simple-minded iteration, with constant temporal step h.
for i=2:n,
    'Starting step', i
    rhssol = y{i-1};
    solinit = bvpinit( rhssol, [-5,5] );
    y{i} = bvp5c(@shocktube, @bc, solinit, opts );
end



end

