How to calculate global stream function field?

Recently, i' m trying to calculate 300hpa global stream function field. So i download the ERA5 data 1.5° X 1.5° velocity u,v field（which size is 121 x 240 x 365) and try to calcute psi.But i' comfused about my result and hope u can give any suggestion. Following is my code and figure.

deltaY = repmat(gsw_distance(double([lon(1) lon(1)]),double([lat(1) lat(2)])),121,1); % meridional 1.5°distance

% cause any 1.5° meridional distance is the same

deltaX = []; %zonal 1.5°distance

for i = 1:length(lat)

deltaX(:,i) = gsw_distance(double(lon),double(repmat(lat(i),length(lon),1)));

end

deltaX = transpose(repmat(deltaX(1,:),240,1));deltaX([1 end],:) = 0;

deltaX = repmat(deltaX,[1 1 365]);

deltaY = repmat(deltaY,[1 1 365]);

index = find(lat == 0);

cy1 = cumsum(u(index-1:-1:1,1,:).*deltaY(index-1:-1:1,:,:));

cy2 = cumsum(u(index:end,1,:).*deltaY(index:end,:,:));

psi1 = cy1(:,ones(size(u,2),1),:) - cumsum(v(index-1:-1:1,:,:).*deltaX(index-1:-1:1,:,:));

psi2 = cy2(:,ones(size(u,2),1),:) - cumsum(v(index:end,:,:).*deltaX(index:end,:,:));

psi = [psi1(end:-1:1,:,:);psi2];

[Dushantha Sandaruwan Jayarathna]

Hi, have you got any sollution or understanding about your work , can I use your given code for calculating psi?