System of partial delay differential equations Ecological models Bifurcation analysis Method of lines Fitted operator finite difference method Convergence analysis a b s t r a c t We consider a system of two coupled partial delay differential equations (PDDEs) describing the dynamics of two cooperative species. The original system is reduced to a system of ordinary delay differential equations (DDEs) obtained by applying the method of lines. Then we construct a fitted operator finite difference method (FOFDM) to solve this resulting system. The model considered in this paper is very sensitive to small changes in the parameters associated in with the model. Depending on the values of these parameters, the solution can be stable, periodic and/or aperiodic. Such behavior of the solution is exploited via the proposed FOFDM. This FOFDM is analyzed for convergence and it is seen that this method is unconditionally stable and has the accuracy of O(k + h 2 ), where k and h denote time and space step-sizes, respectively. Some numerical results confirming theoretical observations are also presented. These results are comparable with those obtained in the literature.