The purpose of the research is to find the numerical solutions to the system of time dependent nonlinear parabolic partial differential equations (PDEs) utilizing the Modified Galerkin Weighted Residual Method (MGWRM) with the help of modified Bernstein polynomials. An approximate solution of the system has been assumed in accordance with the modified Bernstein polynomials. Thereafter, the modified Galerkin method has been applied to the system of nonlinear parabolic PDEs and has transformed the model into a time dependent ordinary differential equations system. Then the system has been converted into the recurrence equations by employing backward difference approximation. However, the iterative calculation is performed by using the Picard Iterative method. A few renowned problems are then solved to test the applicability and efficiency of our proposed scheme. The numerical solutions at different time levels are then displayed numerically in tabular form and graphically by figures. The comparative study is presented along with L2 norm, and L∞ norm.