%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% The first order system of ODE representing the Blasius zeta functions %%%
%%% of the boundary layer along an impulsively started circular cylinder. %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Notations:   x(1) = zeta0,   x(2) = zeta0',   x(3) = zeta0''
%              x(4) = zeta1,   x(5) = zeta1',   x(6) = zeta1''
% The initial condition for zeta0'' at eta = 0 is 2/sqrt(pi) and the one
% for zeta1'' is 2/sqrt(pi)*(1+4/(3*pi)). The other four initial conditions
% are zero.

% Return the derivatives of the x-values
function xp = zetasys(eta, x)
xp = [x(2);
      x(3);
      -2*eta*x(3);
      x(5);
      x(6);
      -2*eta*x(6)+4*x(5)+4*(x(2)^2-x(1)*x(3)-1)];
