Fourth Order Central Differencing in Matlab?

I have been tasked with using a Fourth Order Central Differencing scheme to compute the z-vorticity component (Wz = dv/dx - du/dy) within Matlab; given 3D flow field data for velocities (u,v,w) in a 64^3 point domain plus 4x halo / ghost nodes (2x at the start and 2x at the end) for each dimension.

This problem stems from a 3D Taylor Green Vortex simulation for a cube of size 0<x0,x1,x2<2pi with periodic boundary conditions applied in all directions (time-step of dt=0.005s).

I understand the general concept of central differencing but I do not understand the application to this particular task and the appropriate Matlab syntax so any help would be greatly appreciated.

[J. Alex Lee]

Your question doesn't really indicate to potential helpers where you are stuck, but since I am generally interested in FDM, I can answer at least 2 sub-questions:

1. how to compute coefficients for 4th order differences and
2. efficient matlab code to generate the FD coefficients matrix

For 1., wikipedia is helpful for uniform grids https://en.wikipedia.org/wiki/Finite_difference_coefficient

For non-uniform grids, check out Fornberg, Bengt (1988). "Generation of finite difference formulas on arbitrarily spaced grids". Mathematics of Computation, Volume 51, Number 184, Pages 699-706

For 2., in my experience the most time consuming aspect (for explicit time stepping) is to set up the matrix of FD coefficients, which should be done correctly with sparse(). Here is how I have done for 2D uniform grid for 2nd order in the past (not a full example)

A = N*M; % N x M x P cube
for k = 1:N*M % loop through index
    
    % the subscripts of the central node
    [i,j] = ind2sub([N,M],k);
    
    % for internal nodes, define the cross of nodes centered on (i,j)
    p = sub2ind([N,M],...
	[i,i+1,i-1,i  ,i  ],...
	[j,j  ,j  ,j+1,j-1]);
    
    % make corresponding coefficients according to spatial derivatives of PDE
    s(1) = ...
    s(2) = ...
    s(3) = ...
    s(4) = ...
    s(5) = ...
    
    % subscript index into sparse A
    A = A + sparse(k,p,s,N*M,N*M);
    
end