function bv()



function [x,u,up]=bvshoot(f,xl,xr,lbc,rbc,solver,nint,atol,rtol,u0,up0) 

u0t=u0;
up0t=up0;
options=odeset;
options=odeset(options,'AbsTol',atol);
dnewt=10^(-2)*atol;
       % dnewt=10^(-8);
options=odeset(options,'RelTol',rtol);
        %n=1;
n=nint; 
%solvet=@ode113;   % Default solver
solvet=solver;

xb=linspace(xl,xr,n+1)';
xx=xb(1:n); 

% Build a linear initial guess satisfying the bcs 
    u0t=zeros(n,1);
    up0t=zeros(n,1);
    L=xr-xl;
    B=[ lbc(1) lbc(1)+lbc(2)*L; rbc(1)-rbc(2)*L rbc(1) ];
    z=B\[ lbc(3); rbc(3) ];
    u0t=z(1)*(xx-xl)+z(2)*(xx-xr);
    up0t=(z(1)+z(2))*ones(n+1,1); 

R=zeros(2*n,1);
itmax=100;
it=1;
del=10*dnewt;

while (del > dnewt) & (it <= itmax) 
    R(1)=lbc(1)*up0t(1)+lbc(2)*u0t(1)-lbc(3);
    J=zeros(2*n,2*n);
    J(1,1)=lbc(2);
    J(1,2)=lbc(1);
    for k=1:n
       xspan=[xb(k) xb(k+1)];
       v0=[ u0t(k); up0t(k) ];
       [xsteps v]=solvet(@bvpfun,xspan,v0,options);
       nsteps=length(xsteps);
       if k<n
           R(2*k)=u0t(k+1)-v(nsteps,1);
           R(2*k+1)=up0t(k+1)-v(nsteps,2);
           J(2*k,2*k+1)=1;
           J(2*k+1,2*k+2)=1; 
       else
           R(2*n)=rbc(1)*v(nsteps,2)+rbc(2)*v(nsteps,1)-rbc(3);
       end
       ss=sqrt(eps)*norm(v0);
       v10=v0+[ss; 0];
       [xsteps v1]=solvet(@bvpfun,xspan,v10,options);
       n1=length(xsteps);
       h1=(v1(n1,:)-v(nsteps,:))'/ss;
       v10=v0+[0; ss];
       [xsteps v1]=solvet(@bvpfun,xspan,v10,options);
       n1=length(xsteps);
       h2=(v1(n1,:)-v(nsteps,:))'/ss;
       if k<n
           J(2*k:2*k+1,2*k-1)=-h1;
           J(2*k:2*k+1,2*k)=-h2;
       else
           J(2*n,2*n-1)=rbc(1)*h1(2)+rbc(2)*h1(1);
           J(2*n,2*n)=rbc(1)*h2(2)+rbc(2)*h2(1);
       end
    end
    del=norm(R);
    if (it==1) 
        del0=del;
    else
        if del > 100*del0
            error(' Newton failure ')
        end
    end
    w=J\R;
    for k=1:n
        u0t(k)=u0t(k)-w(2*k-1);
        up0t(k)=up0t(k)-w(2*k); 
    end
    it=it+1;
end

if (it >= itmax)
    warning(sprintf(' No convergence in Newton - final residual=%4.3e',del))
end

% Produce final solution

x=[xl];
u=[u0t(1)];
up=[up0t(1)];

for k=1:n
       xspan=[xb(k) xb(k+1)];
       v0=[ u0t(k); up0t(k) ];
       [xsteps v]=solvet(@bvpfun,xspan,v0,options);
       nsteps=length(xsteps);
       x=[x; xsteps(2:nsteps)];
       u=[u; v(2:nsteps,1)];
       up=[up; v(2:nsteps,2)];
end

%  bvpfun

function yp=bvpfun(x,y)
    yp=zeros(2,1);
    yp(1)=y(2);
    yp(2)=f(y(2),y(1),x);
end

end
   