% Script file: bvPr3
%
% Use continuation to solve problem 3 with eta=100 and eta=1000
%

% We begin with eta=.1

eta=.1;
f=@(up,u,x) eta*u*(1+(u^2-1)*up);
xl=0;
xr=1;
lbc=[0; 1; 1.6];
rbc=[0; 1; -1.7];
[x,u,up]=bvshoot(f,xl,xr,lbc,rbc,@ode45,1,10^(-8));
disp(' eta=0.1 ')

for eta=.2:.4:1
    nint=1;
    f=@(up,u,x) eta*u*(1+(u^2-1)*up);
    xg=linspace(xl,xr,nint+1)';
    u0=pchip(x,u,xg);
    u0p=pchip(x,up,xg);
    [x,u,up]=bvshoot(f,xl,xr,lbc,rbc,@ode45,nint,10^(-8),[],u0(1:nint),u0p(1:nint));
    disp(sprintf(' eta=%2.1f',eta))
end

for eta=2:1:30
    nint=floor(eta/2);
    f=@(up,u,x) eta*u*(1+(u^2-1)*up);
    xg=linspace(xl,xr,nint+1)';
    u0=pchip(x,u,xg);
    u0p=pchip(x,up,xg);
    [x,u,up]=bvshoot(f,xl,xr,lbc,rbc,@ode45,nint,10^(-8),[],u0(1:nint),u0p(1:nint));
    disp(sprintf(' eta=%2.0f',eta))
end

for eta=32:2:80
    nint=eta/2;
    f=@(up,u,x) eta*u*(1+(u^2-1)*up);
    xg=linspace(xl,xr,nint+1)';
    u0=pchip(x,u,xg);
    u0p=pchip(x,up,xg);
    [x,u,up]=bvshoot(f,xl,xr,lbc,rbc,@ode45,nint,10^(-8),[],u0(1:nint),u0p(1:nint));
    disp(sprintf(' eta=%3.0f',eta))
end

for eta=84:4:100
    nint=eta/2;
    f=@(up,u,x) eta*u*(1+(u^2-1)*up);
    xg=linspace(xl,xr,nint+1)';
    u0=pchip(x,u,xg);
    u0p=pchip(x,up,xg);
    [x,u,up]=bvshoot(f,xl,xr,lbc,rbc,@ode45,nint,10^(-8),[],u0(1:nint),u0p(1:nint));
    disp(sprintf(' eta=%3.0f',eta))
end

figure
plot(x,u)
title(' u eta=.01 ')

for eta=150:50:1000
    nint=eta/2;
    f=@(up,u,x) eta*u*(1+(u^2-1)*up);
    xg=linspace(xl,xr,nint+1)';
    u0=pchip(x,u,xg);
    u0p=pchip(x,up,xg);
    [x,u,up]=bvshoot(f,xl,xr,lbc,rbc,@ode45,nint,10^(-8),[],u0(1:nint),u0p(1:nint));
    disp(sprintf(' eta=%3.0f',eta))
end

figure
plot(x,u)
title(' u eta=.001 ')