## 3 Results and Discussion  
Equation (33) was solved numerically in MATLAB R2022a using the bvp5c algorithm, described previously by [PERSON] and [PERSON] [23], with the relative and absolute tolerances set to \(1\times 10^{-12}\). The bvp5c finite difference algorithm uses a collocation formula that provides C1 continuity, which is uniformly fifth-order accurate in the solution interval. The residuals for all calculations were confirmed to be smaller than the tolerance of \(1\times 10^{-12}\). For the purpose of comparison, the analytic solution given by Equation (36) was tabulated for the Reference Case given in Table 1 and compared to the numerical solution calculated via the use of MATLAB. The magnitude of the difference in the values of \(u\) between the two solutions across the channel was less than \(4\times 10^{-17}\). The results obtained for the ranges of flow and domain parameters are given in Table 1.  
The column with the Reference Case heading in Table 1 is used to isolate phenomena and contains the fixed values of parameters when a particular parameter is under investigation. The graphical results for the Reference Case appear in all of the parametric studies to facilitate comparison. For instance, to study the effects of changes in \(\alpha\), its values are varied as listed in the horizontal row labelled \(\alpha\) in Table 1, namely \(\alpha\) takes the range of values \(10^{-1}\), \(10^{-0.5}\), \(1,10^{0.5}\), and \(10^{1}\), while other parameters' fixed values are listed in the reference case column, namely \(\theta=45^{\circ}\), \(\mu_{0}=1\), \(p_{0}=1\), \(n=1\), and \(k=10^{-2}\). In the discussion that follows, the effects of the parameters of Table 1 on viscosity and pressure distributions, velocity profiles, shear stress and viscous drag, and resistance to the flow are discussed. Six parametric studies will be presented in the following order: \(\theta,\ \mu_{0},\ p_{0},\ \alpha,\ n\), and \(k\). The results for each parametric study include (1) velocity profiles, \(u(y)\); (2) viscosity profiles, \(\mu(y)\); (3) Darcy friction profiles, \(-\mu(y)\left(\frac{u(y)}{k}\right)\); (4) the net force on a differential fluid element due to variation in the shear strain rate across the fluid element, \(\mu\left(\frac{d^{2}u}{dy^{2}}\right)\); (5) the net force on a differential fluid element due to variation in the viscosity across the fluid element, \(\left(\frac{du}{dy}\right)\left(\frac{du}{dy}\right)\); and (6) the dimensionless volumetric flowrate in the channel, \(Q\). Items (3)-(5) show the spatial variation in each of the forces present in the dimensionless force balance given by Equation (19) and must sum to \(-\sin(\theta)\), the component of the fluid's weight acting along the incline. Presenting the results in this manner aids in the development of an understanding of the physical phenomena involved in the resulting velocity profiles. Note that in some instances, semi-logarithmic axes are used to better illustrate the results when the quantity in question varies over orders of magnitude. A visual display of the ranges of Table 1 is shown in Figure 2 (below) to demonstrate the increase and decrease in the chosen values for the parameters relative to each other and how they influence the dimensionless volumetric flowrate in the channel, \(Q\). Figure 2 shows that an increase in \(\mu_{0}\) results in a linear decrease in \(Q\). It also illustrates the expected largest rate of increase in \(Q\) with increasing permeability, \(k\).