function bergschrund()
% solve ode for steady glacier thickness with ice flux q = x*(1-x/2)

%% Parameters
n = 3;
mu = 0.1;

%% Numerical solutions of rescaled boundary layer problems
opts = odeset('reltol',1e-6);
xinf = 10; xeps = 1e-4; Heps = 2^(n/(2*n+2))*(n+2)^(1/(2*n+2))*xeps^(1/2); % analytical approximation of singular behaviour
[xh2,Hh2] = ode45(@(x,H) ((n+2)*x.*H.^(-(n+2))).^(1/n)-1 , [xeps xinf] , Heps , opts);
xinf = 10; Hinf = ((n+2)*xinf).^(1/(n+2)); % analytical approximation of far field behaviour
[xh1,Hh1] = ode45(@(x,H) 1-((n+2)*x.*H.^(-(n+2))).^(1/n) , [xinf 0] , Hinf , opts);
Hhstar = Hh1(end),

%% Numerical solution
opts = odeset('reltol',1e-6);
xeps = 1e-4; Heps = mu^(-n/(2*n+2))*2^(n/(2*n+2))*(n+2)^(1/(2*n+2))*xeps^(1/2); % analytical approximation of singular behaviour at x = 2
[xn,Hn] = ode45(@ode_fun,[2-xeps 0],Heps,opts);
function dHdx = ode_fun(x,H)
    dHdx = (1-((n+2)*x.*(1-x/2).*H.^(-(n+2))).^(1/n))/mu;
end

%% Outer solution
x = linspace(0,2,100);
H = ((n+2)*x.*(1-x/2)).^(1/(n+2));

%% Inner solutions
x1 = mu^((n+2)/(n+1))*xh1;
H1 = mu^(1/(n+1))*Hh1;
x2 = 2-mu^((n+2)/(n+1))*xh2;
H2 = mu^(1/(n+1))*Hh2;

%% Plot solutions
figure(1); clf; set(gcf,'Paperpositionmode','auto','units','centimeters','position',[2 4 12 8]);
set(gcf,'DefaultAxesFontSize',12,'DefaultTextFontSize',12); 
    plot(xn,Hn,'b-','linewidth',2);
    hold on;
    plot(x,H,'k--',x1,H1,'r--',x2,H2,'r--');
    xlabel('$x$','interpreter','latex');
    ylabel('$H$','interpreter','latex');
    xlim([0 2]);
    
% %% Save figure    
% print(gcf,'-depsc2',mfilename,'-loose');

%% Plot other numerical solutions of boundary layer at x = 0
figure(2); clf; set(gcf,'Paperpositionmode','auto','units','centimeters','position',[2 4 12 8]);
set(gcf,'DefaultAxesFontSize',12,'DefaultTextFontSize',12); 
    xinf = 10;
    x = linspace(0,xinf,100);
    plot(x,((n+2)*x).^(1/(n+2)),'k--','linewidth',2); % far-field behaviour
    hold on;
    xinf = 10; H0s = [0.9:.02:1.1]*Hhstar; % guesses at initial value
    for H0 = H0s
    [xh1,Hh1] = ode45(@(x,H) 1-((n+2)*x.*H.^(-(n+2))).^(1/n) , [0 xinf] , H0 );
    plot(xh1,Hh1,'-');
    end
    xlabel('$\hat{x}$','interpreter','latex');
    ylabel('$\hat{H}$','interpreter','latex');
   
% %% Save figure  
% print(gcf,'-depsc2',[mfilename,'2'],'-loose');
    
end