I actually solve it. For anyone looking for the same problem:
function [y,dx]=cov_der(x)
%simple example 
%x=dlarray(and(3,3)); 
%[y,dx]=dlfeval(@cov_der,x)
y=zeros(size(x,2),size(x,2)); 
z=size(x,1); 
y=dlarray(y); 
for  i=1:size(x,2)
    for j=1:size(x,2)
    y(i,j)=sum(x(:,i).*x(:,j))./z;
    end
end
dx=zeros(size(x,2)*size(x,2),size(x,1),size(x,2)); 
index=1; 
for i=1:size(x,2)
    for j=1:size(x,2)
    dx(index,:,:)=dlgradient(y(i,j),x,'EnableHigherDerivatives',true);
    index=index+1; 
    end
end
end


