How to compute Jacobian Matrix for coordinate transformation

Hello everyone.

this is my problem I have 4 matrices that represent two different meshes with the same discretization (101x41) X and Y are matrices of a rectangular mesh, i.e they have identical rows. R and Z are matrices of a mesh on a poloidal plane.

I want to compute the Jacobian of the transformation that led from (X,Y) to (R,Z). I've started writing some code but I am getting confused on how to proced because I actually don't now how to compute the derivatives dR/dX dR/dY dZ/dX dZ/dY

I compute the differential on X and Y on a single row and then I copy the result to the other rows since are identical, but then when it comes to compute dR and dZ I don't know hot to proceed. Can anybody help me out?

Here is what I've done so far.

 R=load('RmeshSD2');
 Z=load('ZmeshSD2');
 X=load('XmeshSD2');
 Y=load('YmeshSD2');

 dX=[];
 dR=[];
 dY=[];
 dZ=[];
 IMX=size(Y,1)
 IMY=size(Y,2)
 % return

 %%dX
 for i=2:IMX-1
      dX(i)=(X(1,i)-X(1,i-1));
 end
 dX(1)=interp1(dX(2:IMY-1),X(1,1),'spline','extrap');
 dX(IMX)=interp1(dX(2:IMX-1),X(1,end),'spline','extrap');

 for I=1:IMX
     AX(I,:)=dX;
 end
 dX=AX;

 %%dY
 for i=2:IMY-1
      dY(i)=(Y(1,i)-Y(1,i-1));
 end
 dY(1)=interp1(dY(2:IMY-1),Y(1,1),'spline','extrap');
 dY(IMY)=interp1(dY(2:IMY-1),Y(1,end),'spline','extrap');

 for I=1:IMX
     AY(I,:)=dY;
 end
 dY=AY;

 return

 %%dR
 dR=diff(R);
 %dR(1)=interp1(dR(2:IMX-1),R(1,1),'spline','extrap');
 dR(IMX,:)=interp1(dR(1:IMX-1,:)',R(end,:)','spline','extrap');