function jefferyhamel( Rey, nsteps )
%
% Set up and solve the Jeffery-Hamel BVP by continuation in the Reynolds
% number.
%
    function dy = JHFlow( t, y )
        dy = zeros(size(y));
        dy(1,:) = y(2,:);
        dy(2,:) = y(3,:);
        dy(3,:) = y(4,:);
        dy(4,:) = -2*y(2,:).*y(3,:) - 4*y(3,:);
    end

    function zer = JHBC( ya, yb )
        zer = [ ya(1)
                ya(3)
                yb(2)
                yb(1)+2*Reyl/3 ];
    end

    Reylocal = linspace( 1, Rey, nsteps );

    coarseoptions = bvpset('Vectorized','on', 'RelTol', 1.0e-3, 'AbsTol', 1.0e-3 );

    options = bvpset('Vectorized','on', 'RelTol', 1.0e-7, 'AbsTol', 1.0e-7, 'Nmax', 50000 );

    sol = bvpinit( linspace(0,pi/2,11), [1;0;0;0] );
    
    for i=1:nsteps-1,
        Reyl = Reylocal(i);
        sol = bvp5c( @JHFlow, @JHBC, sol, coarseoptions );
    end
    
    Reyl = Rey;
    sol = bvp5c( @JHFlow, @JHBC, sol, options ); % tight options on last one
        
    t = RefineMesh( sol.x, 20 );
    [y,dy] = deval( sol, t );
    res = dy - JHFlow( t, y );
    
    figure(1),semilogy( t, abs( res(4,:) ), 'k-' ),set(gca,'fontsize',16),xlabel('t','fontsize',16),ylabel('abs(residual)','fontsize',16)
    figure(2),plot( t, y(2,:), 'k-' ), set(gca,'fontsize',16),xlabel('t','fontsize',16), ylabel('f(t) = F''(t)','fontsize',16)


end