%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Blasius second function for impulsively started boundary layer flow. %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% This really returns function zeta_1', which is the scaled second order
% contribution to the velocity profile.  Argument eta may be an array.

% start of the function definition
function z1p = zeta1p(eta)

% constant
c1=1/sqrt(pi);

% store the values of the exponential
ex=exp(-eta.^2);

% store the values of the complementary error function
erc=erfc(eta);

% compute function zeta_1' according to the result of Blasius
z1p=-3*c1*eta.*ex.*erc+.5*(2*eta.^2-1).*erc.^2+2/pi*ex.^2+c1*eta.*ex...
+2*erc-4/(3*pi)*ex+(3+4/(3*pi))*c1*(eta.*ex-.5/c1*(2*eta.^2+1).*erc);
