function pb21_ODE(time_end)
% pb21(time_end) to solve the problem 22.21
% time_end : calculation time to end (0 ~ time_end)
% example 
%  >> pb21(8)
%  The time to height 0 : 6.739 m
[t,y] =ode23(@h,[0 time_end],6,[],0.25,9.81,1);
time=t(1:find(y(:,1)<0));
h=y(1:find(y(:,1)<0));
fprintf('The time to height 0 : %1.3f s\n',time(length(time))); 
plot(time,h)
xlabel('Time [s]'); ylabel('Height [m]')
grid on
end

function dh=h(t,y,d,g,e)
h=[6 5 4 3 2 1 0];
A=[1.17 0.97 0.67 0.45 0.32 0.18 0];
Ah=interp1(h,A,y);
dh=[-pi*d^2*sqrt(2*g*(y+e))/(4*Ah)];
end