The study introduces a computational method that combines Legendre polynomials with Gauss-Lobatto points to solve nonlinear coupled differential equations, focusing on the Williamson fluid model with the existence of mixed convection and permeability under mixed boundary conditions. The nonlinear governing equations were transformed to ordinary differential equation (ODE) from partial differential equation (PDE), applying the appropriate similarity conversions. By using Legendre polynomials as trial functions and collocating residual equations with Gauss-Lobatto points, the system is solved with Mathematical software. The technique was validated by comparing the obtained solution with an existing literature and further validation was done with Runge-Kutta of order 4 via shooting method. Validation against the Shooting Runge-Kutta method showed minimal discrepancies, confirming the method’s accuracy. Graphical analysis indicated that an increase in the Grashof number enhances velocity, while higher porosity raises temperature but reduces fluid velocity. This approach offers an efficient and precise solution for complex nonlinear equations, with broader potential applications in fluid dynamics.