In this paper, we have used a streamfunction-vorticity (ψ − ξ) formulation to investigate the problem of 2-D unsteady viscous incompressible flow in a driven square cavity with moving top and bottom walls. We used this formulation to solve the governing equations along with no-slip and slip wall boundary conditions. A general algorithm was used for this formulation in order to compute the numerical solutions for the flow variables: streamfunction ψ, vorticityfunction ξ for low Reynolds numbers Re ≤ 50. We have executed this with the aid of a computer programme developed and run in C++ compiler. We have proved the stability and convergence of the numerical scheme using matrix method. Following this, the stability conditions obtained for the time and space steps have been used in numerical computations to arrive at the numerical solutions with desired accuracy. Also investigated was the variation of u and v components of velocity at points closed to left side, top and bottom walls of the square cavity at different time levels for Reynolds number 50. Based on the numerical computations of u-velocity, we found: (i) the u-velocity is almost same near the top and bottom wall of the square cavity above and below the geometric center for Re=15 and 30; (ii) the absolute value of the u-velocity first decreases, then increases, and finally decreases, near top and bottom wall of the square cavity. Based on the numerical computations of the v-velocity, we found: (i) the absolute value of the v-velocity increases uniformly as time level increases' (ii) the v-velocity at the right wall of the square cavity attains a maximum value, and then, decreases towards the geometric center. The numerical solutions of the vorticity vector ξ at Re=50 for different times along the horizontal line through geometric center of the square cavity indicate the following: (i) the vorticity contours decreased in between the left wall boundary to the midpoint of the square domain; (ii) it, then, increased in between the midpoint to the right wall boundary. Based on the numerical solutions of streamfunction ψ, we found: (i) the size of streamline contours decreased when the time increased; (ii) two streamline contours are generated, one above the geometric center, and the other, below the geometric center in clockwise direction.