I am trying to solve the Falkner Skan equation for wedge flows numerically, by first reducing it to a second order system and then solving a tridiagonal matrix equation iteratively till solution converges to a required tolerance. Essentially we have: $$ f''' + ff'' + \beta(1-f'^2) = 0\\ f(0)= f'(0)=0\\ f'(\infty)=1\\ $$ Here, we take $$ f' = G$$ And the equation becomes $$ G'' + G' * \int_{0}^{\eta} G d\eta + (1-G^2)\beta = 0\\ G(0) = f'(0)= 0\\ G(\infty=n) = 1$$ The nonlinear part in the RHS has been approximated as: $$(G^2 - 1) = 2G_{{j}_{old}}G_{j} - (G_{{j}_{old}})^2 - 1$$ G_{j_{old}} is the previous guess in the iterative scheme. We then use the central differences rules for G and write the equation in terms of $G_{j-1}, G_j, G_{j+1}$ over a non-uniform grid of n grid points from 0 to $\eta_n$ i.e. $\eta_j$ where j varies from 0 to n. The equation is something like this, $$a_j G_{j-1} + b_j G_j +c_j G_{j+1} = -(G_{j_{old}}^2+1)\beta$$ And the tridiagonal matrix equation looks like this, \begin{equation} \begin{bmatrix} b_2 & c_2 & 0 & 0 & - & -\\ a_3 & b_3 & c_3 & 0 & - & - \\ 0 & a_4 & b_4 & c_4 & 0 & -\\ 0 & 0 & a_5 & b_5 & c_5 & -\\ 0 & 0 & 0 & - & - & - \\ - & - &- & a_{n-2} & b_{n-2} & c_{n-2}\\ - & - & - & - & a_{n-1} & b_{n-1} \end{bmatrix} \times \begin{bmatrix} G_2 \\ G_3 \\ G_4\\ G_5\\-\\-\\ G_{n-1} \end{bmatrix} = \begin{bmatrix} -\beta(1+G^{old^2}_{2})-a_2 G_1 \\ -\beta(1+G^{old^2}_{3}) \\ -\beta(1+G^{old^2}_{4})\\ -\beta(1+G^{old^2}_{5})\\-\\-\\ -\beta(1+G^{old^2}_{n-1}) - c_{n-1} G_n \end{bmatrix} \end{equation} And, now we use the Thomas algorithm to invert the following equation and find $G (or f')$ from this, $$ AG = R $$
I solve this in Fortran and my results look like this:
They are supposed to look like this
I am unable to understand why this works for positive values of $\beta$ but fails for negative and zero.