function [ y ] = burgerstmol( tf )
%BURGERSTMOL Solve Burgers' equation by TMOL
%   Use bvp4c and Backward Euler
%   y = burgerstmol( tf );  v = u_x
%   u_t = -(u^2/2)_x + nu* v_x

nu = 2/100;
% Initial and boundary values defined below
    function y = inivals( x )
        k = find(x<=0.1);
        m = find(x>0.1);
        y = [ones(size(x(k))),zeros(size(x(m)))
             zeros(size(x))];
    end
    function res = bc( ya, yb )
        res = [ya(1)-1
               yb(1)];
    end

inimesh = [0,0.09,0.099,0.101,0.11, 1., 2, 3, 4];
solinit = bvpinit( inimesh, @inivals );

n = 80;
h = tf/n; % Apathetic choice.
y = cell(n);

% Initial values used on the first step.
    function yp = advectionbvp1( x, y )
        rhsx = inivals(x);
        yp = zeros(size(y));
        yp(1,:) = y(2,:);
        yp(2,:) = y(1,:).*y(2,:)/nu + (y(1,:) - rhsx(1,:))/h/nu;
    end

% Take the first step from the initial condition
solinit = bvpinit( inimesh, @inivals );
y{1} = bvp4c( @advectionbvp1, @bc, solinit );

% Subsequent steps use the previous solution
% as initial conditions, and also as the right-hand side.
% This is passed to the following function by nested scoping.
    function yp = advectionbvp( x, y )
        [yo,yox] = deval( rhssol, x );
        rhsx = yo(1,:); 
        yp = zeros(size(y));
        yp(1,:) = y(2,:);
        yp(2,:) = y(1,:).*y(2,:)/nu + (y(1,:) - rhsx(1,:))/h/nu;
    end

% Simple-minded iteration, with constant temporal step h.
for i=2:n,
    rhssol = y{i-1};
    solinit = bvpinit( rhssol, [0,4] );
    y{i} = bvp4c(@advectionbvp, @bc, solinit );
end

end

