problem implementing central finite difference on multicore computer?

Hi,

I am trying to implement the central finite difference scheme for use on a multicore computer. The computer has 12 cores.

Currently I have a code that takes several days to run and any help in improving the performance would be helpful.

My finite difference function looks like this:

 %nx,ny,nz = number of grid points in x,y,z direction
 %dx = grid spacing
 %v = velocity model
 %ricker = ricker wavelet
 %nt = number of time samples
 %dt = number of time samples
 %xs,ys,zs = source position
% receivers = receiver positions in y direction

 function u = fd3d(nx,ny,nz,dx,v,ricker,nt,dt,xs,ys,zs,receivers)

 p0 = zeros(nx,ny,nz); p1 = p0; p2 = p0;
 cns = 6*(dt*v).^2;

 u = zeros(length(receivers),nt);
 alpha=dt/dx;

 for it = 1:nt
    p2 = 2*p1 - p0 + cns.*del2(p1,dx);
    p2(xs,ys,zs) = p2(xs,ys,zs) + ricker(it);

    %ABC BOUNDARY CONDITIONS

    p2(1,:,:) = p1(1,:,:) + alpha.*v(1,:,:).*(p1(2,:,:)-p1(1,:,:));

    p2(end,:,:) = p1(end,:,:) + alpha.*v(end,:,:).*(p1(end-1,:,:)-p1(end,:,:));

    p2(:,1,:) = p1(:,1,:) + alpha.*v(:,1,:).*(p1(:,2,:)-p1(:,1,:));

    p2(:,end,:) = p1(:,end,:) + alpha.*v(:,end,:).*(p1(:,end-1,:)-p1(:,end,:));

    p2(:,:,1) = p1(:,:,1) + alpha.*v(:,:,1).*(p1(:,:,2)-p1(:,:,1));

    p2(:,:,end) = p1(:,:,end) + alpha.*v(:,:,end).*(p1(:,:,end-1)-p1(:,:,end));

    a = p2(xs,receivers,zs); extract seismogram from current timestep
    u(:,it) = a';  %store seismogram 
    p0 = p1; p1 = p2; %update pressure fields

 end

The function that takes the longest time to compute is the builtin del2 function. I was hoping if anyone had any suggestion on a function that can be used instead of del2 that works together with the parallel computation toolbox.