close all; clear all; clc

% calculate h(x) using ODE45 integrator in Matlab
% h = h1, dh/dx=h2
% units: mm, 

initvals(1)=input('h(x=0):');
initvals(2)=input('dh/dx(x=0):');
xlo=0;
xmax=10;
xstep=0.5;
xlim=0.;

while xlim <= xmax
    xlim=xlim+xstep;
    [x,h]=ode45(@DEs,[xlo xlim],initvals);
    psol=pi/2*exp(-x/2);
    plot(x,h(:,1),x,psol);
    last=length(x);
    slope(1)=(h(last)-h(last-1))/(x(last)-x(last-1));
    slope(2)=(h(last-1)-h(last-2))/(x(last-1)-x(last-2));
    if (abs(slope(1)) > abs(slope(2))) && (xlim < xmax)
        initvals(1)=input('new value of h(x=0):');
        xlim=xstep;
    end
end

%save results once sufficiently converged
savefile='results.txt';    
fileID=fopen(savefile,'w');
hh=h(:,1);
for i=1:length(x)
    fprintf(fileID,'%d %d\n',x(i),hh(i));
end
fclose(fileID);
disp('results saved to results.txt')

function [dhdx]=DEs(x,h)
dhdx=zeros(2,1);
lc=2; % cappilary length
dhdx(1)=h(2);
dhdx(2)=(1/lc^2)*h(1)*(1+h(2)^2)^(3/2);
end

