function sol = shocknozzle4c
%
% Shock wave in nozzle flow
% Example 1.17 in Ascher, Mattheij, and Russell
%
% Attempted solution by bvp4c---which fails.
% Replace bvp4c with bvp5c in what follows, and it apparently succeeds.
% RMC July 2011
%
    function du = shocknozzle(x,u)
        A = 1 + x.*x;
        dA = 2*x;
        du = zeros(size(u));
        du(1,:) = u(2,:);
        du(2,:) = ((1+gamma/2-epsilon*dA).*u(1,:).*u(2,:) ...
            - u(2,:)./u(1,:) - (dA./A).*(1-(gamma-1)*u(1,:).^2/2 )) ...
            /epsilon ./A ./u(1,:);
    end

    function y = shockbc( ua, ub )
        y = [ua(1)-0.9129
             ub(1)-0.375];
    end

solinit = bvpinit( linspace(0,1,10), [0.9;0] );
        
gamma = 1.4;
% simple continuation in epsilon
n = 40;
pow = linspace(1,8,n); % Chebfun with 512 died at 151, or epsilon = 4.224e-3
maxres = zeros(n,1);
xi = linspace(0,1,1001);
for i=1:n,
    epsilon = 4.792 * 10^(-pow(i));
    disp( [i,epsilon] );
    opts = bvpset('reltol',min(1.0e-6,epsilon/5),'NMax',10000);
    sol = bvp5c(@shocknozzle,@shockbc,solinit,opts);
    solinit = bvpinit( sol, [0, 1] );
    [u,du] = deval( sol, xi );
    res = du - shocknozzle(xi,u);
    maxres(i) = max(max(abs(res)));
end
% Polish the last solution a bit
opts = bvpset('reltol',min(1.0e-6,epsilon/10),'NMax',50000);
sol = bvp5c(@shocknozzle,@shockbc,solinit,opts);
[u,du] = deval( sol, xi );
res = du - shocknozzle(xi,u);
figure(1),plot(xi,u(1,:),'k');
xlabel('x','fontsize',16);
ylabel('solution','fontsize',16);
set(gca,'fontsize',16);
figure(2),semilogy( abs(res(2,:)), 'k' );
xlabel('x','fontsize',16);
ylabel('residual');
set(gca,'fontsize',16);
figure(3),semilogy( sol.x(2:end), diff(sol.x), 'k.' );
set(gca,'fontsize',16);
xlabel('x','fontsize',16);
ylabel('mesh widths','fontsize',16)
end