This paper proposes a method to numerically study viscous incompressible two-dimensional steady flow in a driven square cavity with heat and concentration sources placed on its side wall. The method proposed here is based on streamfunction-vorticity (Ψ-ξ) formulation. We have modified this formulation in such a way that it suits to solve the continuity, x and y-momentum, energy and mass transfer equations which are the governing equations of the problem under investigation in this study. No-slip and slip wall boundary conditions for velocity, temperature and concentration are defined on walls of a driven square cavity. In order to numerically compute the streamfunction Ψ, vorticityfunction ξ , temperature θ, concentration C and pressure P at different low, moderate and high Reynolds numbers, a general algorithm was proposed. The sequence of steps involved in this general algorithm are executed in a computer code, developed and run in a C compiler. We propose that, with the help of this code, one can easily compute the numerical solutions of the flow variables such as velocity, pressure, temperature, concentration, streamfunction, vorticityfunction and thereby depict and analyze streamlines, vortex lines, isotherms and isobars, in the driven square cavity for low, moderate and high Reynolds numbers. We have chosen suitable Prandtl and Schmidt numbers that enables us to define the average Nusselt and Sherwood numbers to study the heat ad mass transfer rates from the left wall of the cavity. The stability criterion of the numerical method used for solving the Poisson, vorticity transportation, energy and mass transfer has been given. Based on this criterion, we ought to choose appropriate time and space steps in numerical computations and thereby, we may obtain the desired accurate numerical solutions. The nature of the steady state solutions of the flow variables along the horizontal and vertical lines through the geometric center of the square cavity has been discussed and analyzed. To check the validity of the computer code used and corresponding numerical solutions of the flow variables obtained from this study, we have to compare these with established steady state solutions existing in the literature and they have to be found in good agreement.