nuair = 1.6e-5; % kinematic viscosity m^2/s (same unit as alpha)
nuwater = 1e-6;

if(Substance == 1)
    nu = nuair; 
elseif(Substance == 2)
    nu = nuwater;
end

% pressure gradient (if convection is used)
dpdx = 10; 
% A = 0 -> pure conduction problem (heat eqn)
% A = 1 -> conduction and convection problem (CD eqn)
A = 1;
% velocity
hw = Ly/2; % channel half width
U = A*(dpdx/(2*rho*nu))*(hw^2)*(1-((Y-hw)./hw).^2); % parabolic U=U(y)
V = 0*U;