function figureS1
figpref(4)

%-- Load SE data and extract important parameters

a = load('newdata/fmode_3_24_T0_new.mat');   % best agreement is with hT0 data
f = a.f; x = a.x; n = a.n;
T = a.tension; h = a.thickness; L = a.length; W = a.width; Delta = a.delta;
Y = 1; nu = 0.4; E = Y/h; B = E*h^3/(12*(1-nu^2));

N0 = 3; N1 = 24;
nbar = @(x) 1./((1/N0) + (1/N1-1/N0)*(x-0)/(L-0));
lambdabar = @(nbar) W./nbar;
Ksub = @(x) B*(2*pi*nbar(x)/W).^4;

n0 = 5:20;
x0 = zeros(size(n0)); x1=x0; xm=x0;
cols = parula(length(n0)+4);

OPTION = 1; %gives best agreement

if OPTION == 1
    % option 1: choose intervals [n0,(n0+n1)/2] = [n0,n0+1/2];
    n1 = n0+1/2; nm = n0-1/2;
    options.Display = 'off';
    for i = 1:length(n0)
        x0(i) = fsolve(@(x) nbar(x)-n0(i),0,options);
        x1(i) = fsolve(@(x) nbar(x)-n1(i),0,options);
        xm(i) = fsolve(@(x) nbar(x)-nm(i),0,options);
    end
end
if OPTION==2
    % option 2: choose intervals [lambda0, (lambda0+lambda1)/2
    options.Display = 'off';
    lambda0 = lambdabar(n0); 
    lambda1 = 0.5*(lambdabar(n0)+lambdabar(n0+1));
    lambdam = 0.5*(lambdabar(n0)+lambdabar(n0-1));
    n1 = fsolve(@(n) lambdabar(n)-lambda1,n0,options);
    nm = fsolve(@(n) lambdabar(n)-lambdam,n0,options);
    for i = 1:length(n0)
       x0(i) = fsolve(@(x) lambdabar(nbar(x))-lambda0(i),0,options);
       x1(i) = fsolve(@(x) lambdabar(nbar(x))-lambda1(i),0,options);
       xm(i) = fsolve(@(x) lambdabar(nbar(x))-lambdam(i),0,options);
    end
end
if OPTION==3
    % option 3: choose actual mid-point
    options.Display = 'off';
    for i = 1:length(n0)
        x0(i) = fsolve(@(x) nbar(x)-n0(i),0,options);
        x1(i) = 0.5*(x0(i) + fsolve(@(x) nbar(x)-n0(i)-1,0,options));
        xm(i) = 0.5*(x0(i) + fsolve(@(x) nbar(x)-n0(i)+1,0,options));
    end
    n1 = nbar(x1); nm = nbar(xm);
end

Kstar = Ksub(x0);
delta = 2*pi*(B./(Kstar*W^4)).^(1/4);
Ls = delta.*(x1-x0).*sqrt(Kstar/T);
Ll = delta.*(xm-x0).*sqrt(Kstar/T);

%--

% Build initial guess

kappa = @(X) 10*X;
solinit.x = linspace(0,0.5/10);
solinit.y = [linspace(0.1,0.2);linspace(1,1);linspace(0.2,0.1);linspace(1,1);...
              linspace(0.1,0.2);linspace(1,1)];
sol = bvp4c(@(X,y) odefun(X,y,@(X) kappa(X)), @(ya,yb) bcfun(ya,yb), solinit);

alpha = logspace(1,log10(0.5),20); %I need the first ones to get a good initial guess
for i = 2:length(alpha)
     kappa = @(X) alpha(i)*X;
     solinit.x = alpha(i-1)*sol.x/alpha(i); 
     solinit.y = sol.y;
     sol = bvp4c(@(X,y) odefun(X,y,@(X) kappa(X)), @(ya,yb) bcfun(ya,yb), solinit);
end
X0 = 2*alpha(end)*sol.x; Y0 = sol.y;

% Solve for experimental values
for i = 1:length(n0)
    
    % n0 to n0+1/2 (or lambda0 to lambda0+...)
    kappa = @(X) 0.5*(X/Ls(i));
    solinit.x = X0*Ls(i); solinit.y = Y0;
    sol = bvp4c(@(X,y) odefun(X,y,@(X) kappa(X)), @(ya,yb) bcfun(ya,yb), solinit);
    X = sol.x'; y = sol.y'; 
    t0 = y(:,1); t1 = y(:,2); t2 = y(:,3); 
    
    figure(1)
    plot(x0(i)+sqrt(T/Kstar(i))*X/delta(i),sin(t0).*cos(t1),'color',cols(i+1,:));
    hold on;
    plot(x0(i)+sqrt(T/Kstar(i))*X/delta(i),sin(t0).*sin(t1).*sin(t2),'color',cols(i+2,:));
    plot(x0(i)+sqrt(T/Kstar(i))*X/delta(i),sin(t0).*sin(t1).*cos(t2),'color',cols(i+3,:));
    plot(x0(i)+sqrt(T/Kstar(i))*X/delta(i),cos(t0),'color',cols(i+4,:));

    %n0 to n0-1/2 (or lambda0 to lambda0-...)
    kappa = @(X) -0.5*(X/Ll(i));

    solinit.x = X0(end:-1:1)*Ll(i); solinit.y = Y0;
    sol = bvp4c(@(X,y) odefunm(X,y,@(X) kappa(X)), @(ya,yb) bcfunm(ya,yb), solinit);
    X = sol.x'; y = sol.y'; 
    t0 = y(:,1); t1 = y(:,2); t2 = y(:,3); 
    
    plot(x0(i)+sqrt(T/Kstar(i))*X/delta(i),sin(t0).*cos(t1),'color',cols(i+3,:));
    hold on;
    plot(x0(i)+sqrt(T/Kstar(i))*X/delta(i),sin(t0).*sin(t1).*sin(t2),'color',cols(i+2,:));
    plot(x0(i)+sqrt(T/Kstar(i))*X/delta(i),sin(t0).*sin(t1).*cos(t2),'color',cols(i+1,:));
    plot(x0(i)+sqrt(T/Kstar(i))*X/delta(i),cos(t0),'color',cols(i,:));
    hold on;    
    
    % add SE data to plot
    m = n0(i)-2;
    plot(x{m},f{m}*double(n(m))/(sqrt(-Delta*W^2)/pi),'--','color',cols(i+2,:));

end

hold off;
xlim([40,78]);
xlabel('$x$','interpreter','latex');
ylabel('$k_i A_i$','interpreter','latex');

hold on; 
for i = 1:length(n0)
    plot(x1(i)*ones(100,1),linspace(0,1),'k--','LineWidth',1);
end
hold off;
end


function dy = odefun(X,y,kappa)

t0 = y(1); t1 = y(2); t2 = y(3); t0p = y(4); t1p = y(5); t2p = y(6);
K = kappa(X);

k0 = 2; k1 = -1; k2 = +1; k3 = 0;
F0 = 4*(k0-K)^2;
F1 = 4*(k1-K)^2; 
F2 = 4*(k2-K)^2;
F3 = 4*(k3-K)^2;

dy = [t0p; t1p; t2p; ...
      sin(t0)*cos(t0)*(t1p^2 + sin(t1)^2*t2p^2 +...
         (F3*sin(t2)^2+F2*cos(t2)^2)*sin(t1)^2 + F1*cos(t1)^2 - F0);...
       sin(t1)*cos(t1)*(t2p^2 + F3*sin(t2)^2 + F2*cos(t2)^2-F1)...
            - 2*t0p*t1p*cos(t0)/sin(t0); ...
       sin(t2)*cos(t2)*(F3-F2) - 2*t0p*t2p*cos(t0)/sin(t0) - 2*t1p*t2p*cos(t1)/sin(t1)];
%pause
end

function res = bcfun(ya,yb)

t0L = ya(1); t1L = ya(2); t2L = ya(3); t0pL = ya(4); t1pL = ya(5); t2pL = ya(6);
t0R = yb(1); t1R = yb(2); t2R = yb(3); t0pR = yb(4); t1pR = yb(5); t2pR = yb(6);

res = [t0L-pi/2; cos(t1L)-sin(t1L)*cos(t2L); t1pL*sin(t2L) + t2pL;...
    t2R-pi/4; cos(t0R)-sin(t0R)*cos(t1R); t1pR+t0pR*sin(t1R)];

end

function dy = odefunm(X,y,kappa)

t0 = y(1); t1 = y(2); t2 = y(3); t0p = y(4); t1p = y(5); t2p = y(6);
K = kappa(X);

k0 = -2; k1 = +1; k2 = -1; k3 = 0;
F0 = 4*(k0-K)^2;
F1 = 4*(k1-K)^2; 
F2 = 4*(k2-K)^2;
F3 = 4*(k3-K)^2;

dy = [t0p; t1p; t2p; ...
      sin(t0)*cos(t0)*(t1p^2 + sin(t1)^2*t2p^2 +...
         (F3*sin(t2)^2+F2*cos(t2)^2)*sin(t1)^2 + F1*cos(t1)^2 - F0);...
       sin(t1)*cos(t1)*(t2p^2 + F3*sin(t2)^2 + F2*cos(t2)^2-F1)...
            - 2*t0p*t1p*cos(t0)/sin(t0); ...
       sin(t2)*cos(t2)*(F3-F2) - 2*t0p*t2p*cos(t0)/sin(t0) - 2*t1p*t2p*cos(t1)/sin(t1)];
%pause
end

function res = bcfunm(ya,yb)

t0L = ya(1); t1L = ya(2); t2L = ya(3); t0pL = ya(4); t1pL = ya(5); t2pL = ya(6);
t0R = yb(1); t1R = yb(2); t2R = yb(3); t0pR = yb(4); t1pR = yb(5); t2pR = yb(6);

res = [t0R-pi/2; cos(t1R)-sin(t1R)*cos(t2R); t1pR*sin(t2R) + t2pR;...
    t2L-pi/4; cos(t0L)-sin(t0L)*cos(t1L); t1pL+t0pL*sin(t1L)];

end