The established method of splitting the equations of flow into separate convection and diffusion equations, the solutions of which are superposed to provide an approximate solution, is here presented with a new vorticity diffusion treatment. Test cases which show the accuracy of the method in several simple flows are presented. The method is applied to low Reynolds number flows about circular cylinders with both steady and oscillating wakes. The diffusion method involves a direct exact solution for the diffusion of vorticity at each time step followed by rediscretisation onto a mesh of grid points. Comparison with a finite element solution and with experimental results shows good agreement up to a Reynolds number of 200, the extent of the range at present studied. Two meshes of grid points are used, firstly an exponentially expanding polar mesh extending from the circular body surface over the whole flow field; this is used to obtain the stream function solution of the Poisson's equation for convection by the Cloud-in-Cell method. Secondly, a mesh appropriate to the diffusion solution is used; this is a cartesian mesh of fixed spacing in the wake region outside the boundary part of the polar mesh. At each diffusion time step discrete point vortices are introduced on the boundary of the body to represent a vortex sheet between the no-slip surface and the velocity calculated at the surface at that time. Interaction of the diffusing vortices with the surface of the cylinder is modelled by the reflection of the diffusing vorticity from the surface. In the diffusion process, at the end of each diffusion time step, vorticity is calculated only at those mesh points within the area where the vorticity becomes significant. In the diffusion from the individual vortices the circulation is conserved. The contributions at each mesh point from all the neighbouring vortices are combined and new zero-age vortices are created, provided that the vorticity is above a small cut-off value. The choice of this value has a significant effect on programme efficiency; various values were tried to ensure sufficient accuracy. In the symmetrical cylinder wake studies, the effects of variations of all the numerical parameters were scrutinised. The effects on the forces, surface vorticity and wake geometry were small. These tests also showed that the diffusion time step should be chosen so that diffusion spreads to cover more than two mesh lengths. In the asymmetic model, only the effects of reducing the time step and vorticity cut-off have been checked and found to be insignificant. In the case of steady flow at Reynolds numbers 20 and 40, favourable comparison is made with the results of the finite difference method of Braza et al. (1986) for length of separation bubble, angle of separation, drag coefficient, surface pressure coefficient and surface vorticity. Comparison with experimentally determined quantities also produces good agreement. The oscillating wake was initiated at Reynolds numbers of 100 and 200 by placing asymmetrically in the wake one concentrated vortex. Steady state oscillations were observed after a very short time-span and followed for 3 complete cycles during which the amplitudes of oscillating lift were constant. Good agreement with flow visualisation and with the results of Braza et al. were obtained with the possible exception of the oscillation of the separation angle.