Pee computation

r=0.004;
m = (4*pi*r^3 / 3 * 1000);
A=pi*r^2;
cd=0.5;
rhoair=1.29;
a = 1/2*rhoair*cd*A/m;
b=9.81;

t=0:0.01:5;

y=sqrt(b/a)*tanh(sqrt(a*b)*t);
p=cumtrapz(t,y);

ax=plotyy(t,y,t,p);
hold on
plot(t, sqrt(b/a)*ones(1, numel(t)), 'r:')
hold off
xlabel('time')
legend('velocity [m/s]', 'terminal velocity', 'distance [m]', 'Location', 'SouthEast');
ylabel(ax(1), 'velocity [m/s]');
ylabel(ax(2), 'distance [m]');