function arc_length

N   = 1000;
x   = linspace(0,1,N);
dx  = x(2)-x(1);

% Original path
%------------------------------------------------------------
y   = 1 - x;
plot(x,y); hold on;
for k = 1:N-1
    dy(k) = (y(k+1)-y(k))/dx; 
end

% Modified path
%------------------------------------------------------------
eps = 0.3;
del = eps*x.*(1-x);
y2  = y + del;
plot(x,y2,'r');
for k = 1:N-1
    dy2(k) = (y2(k+1)-y2(k))/dx; 
end

% Compute arc lengths by direct computation
%------------------------------------------------------------
sum = 0;
for k = 1:N-1
   sum = sum + sqrt(1+dy(k)^2)*dx; 
end
arc1 = sum

sum = 0;
for k = 1:N-1
   sum = sum + sqrt(1+dy2(k)^2)*dx; 
end
arc2 = sum

diff_direct = arc2 - arc1

% Compute difference in arc lengths by integral
%------------------------------------------------------------
for k = 1:N-1
    integrand(k) = eps*dy(k)*dy2(k)/sqrt(1+dy(k)^2);
end


sum = 0;
for k = 1:N-1
   sum = sum + integrand(k)*dx; 
end
diff_indirect = sum
