Solving systems of convection diffusion equations

I'm attempting to use MATLAB to solve a system of 2D convection diffusion equations:

dx/dt = Mx + D\nabla^2 x

Where x is a vector and M and D are matrices (I'm likely only trying to solve two equations at once). I'm not entirely sure how to start, or whether MATLAB is well-equipped for this problem. I've tried to use solvepde as follows:

function results = plotgraph()

model = createpde(2);
d1 = [2,2,2,2;0,1,1,0;1,1,0,0;0,0,10,10; 0, 10, 10, 0; 1, 1, 1, 1; 0 0 0 0];
geometryFromEdges(model, d1);
generateMesh(model, 'Hmax', 0.25);
applyBoundaryCondition(model, 'neumann', 'Edge', [1:4], 'g', [0:0]);
specifyCoefficients(model, 'm', 0, 'd', 1, 'c', [1;5], 'a', [-4.385;-5.385;7.306;10.959], 'f', [0;0]);
u0 = @u0fun;
setInitialConditions(model, u0);
results=solvepde(model, 0:5000);

return

Where u0fun defines some initial conditions. However, I quickly get

Warning: Failure at t=6.255833e+02.  Unable to meet integration tolerances without reducing the step size below the smallest value
allowed (1.818989e-12) at time t.

The system also does not evolve how I'd expect it to in short time frames, which makes me suspect I might be doing something wrong. Any suggestions?

[Sanket Karnik]

I understand that you are getting 'Unable to meet integration tolerances without reducing the step size below the smallest value' error when you try to execute above code. Can you please post your 'u0fun' code, so that we can try to use same initial conditions?

Thanks!