function result = K_tau_12(x,y,t)
    m=100;
    [Wweight, Wt] = ClenshawCurtis(0, 100, m);
    [Vweight, Vs] = ClenshawCurtis(0, 100, m);
     
    [v, w] = ndgrid(Vs, Wt); 
    [vweight, wweight] = ndgrid(Vweight, Wweight);
     
    integrand_vals = integrand(x + v, w + y).*exp(t.*(w-v));
    remainder1=exp(t^3/3-t*y)*Vweight*(airy(x+Vs).*Vs.*exp(-t.*Vs));
    remainder2=exp(-t^3/3+t*x)*Vweight*(airy(y+Vs).*Vs.*exp(t.*Vs));
    remainder3 = exp(t*(x-y))*(t^(2)-min(x,y));
    weights = vweight .* wweight; 

    result = real(-sum(integrand_vals .* weights, 'all')+remainder1+remainder2-remainder3);
    %result=real(-sum(integrand_vals .* weights, 'all'));
    %result = real(result*conjugation(x,y,t));
end

function s = integrand(x, y) 
    airy_x = airy(x);
    airy_y = airy(y);
    airy1_x = airy(1, x);
    airy1_y = airy(1, y);
     
    diff = x - y;
    is_equal = abs(diff) ==0;  
     
    s = zeros(size(x));
     
    s(~is_equal) = (airy_x(~is_equal) .* airy1_y(~is_equal) - ...
                    airy1_x(~is_equal) .* airy_y(~is_equal)) ./ diff(~is_equal);
    
     
    s(is_equal) = -x(is_equal) .* airy_x(is_equal).^2 + airy1_x(is_equal).^2;
end

function s = conjugation(x,y,t)
    s= exp(-t*(x-y));
end