Bilinear interpolation in a 2d model

Hello all,

I am writing a particle tracking model on matlab to study the larvae dispersal of mussel. The code works well. However, the particles are moving in 2 dimensions and i used a simple linear interpolation in my code. I would like to change it into a bilinear interpolation. Can someone help me please? find below the part of the code which does the linear interpolation:

% linear interpolation of velocity:
    clear I J I2 J2 I3 J3 Itemp Jtemp I2temp J2temp incx incy u v u1 u2 v1 v2
    
    % find the boxes sup and inf where the particle is:
    I = max(find(xq<x(it,ip)));
    J = max(find(yq<y(it,ip)));
    I2 = find(xq>x(it,ip),1);
    J2 = find(yq>y(it,ip),1);
    
    % Value of the limit boxes in UTM:
    Itemp = xmin + (dx*(I-1));
    Jtemp = ymin + (dy*(J-1));
    I2temp = xmin + (dx*(I2-1));
    J2temp = ymin + (dy*(J2-1));
    
    % index creation
    incx = (x(it,ip)-Itemp)/(I2temp-Itemp);
    incy = (y(it,ip)-Jtemp)/(J2temp-Jtemp);
    
    u1 = MASK(J,I).*U(J,I);
    u2 = MASK(J,I2).*U(J,I2);
    v1 = MASK(J,I).*V(J,I);
    v2 = MASK(J2,I).*V(J2,I);
    
    if MASK(J,I) == 1 & MASK(J,I2) == 1;
        u = u1 + incx*(u2-u1);
    elseif MASK(J,I)== 1 & MASK(J,I2 )== 0;
        u = u1;
    else MASK(J,I)== 0 & MASK(J,I2 )== 1;
        u = u2;
    end
    
    if MASK(J,I) == 1 & MASK(J2,I) == 1;
        v = v1 + incy*(v2-v1);
    elseif MASK(J,I) == 1 & MASK(J2,I) == 0;
        v = v1;
    else MASK(J,I) == 0 & MASK(J2,I) == 1;
        v = v2;
    end