clear all; close all

d=1;            %breite der Quellverteilung
geschw=1;       %Geschwindigkeit im unendlichen

%%%%%%%%%%%%%%%%%
x1=[-d -0.1];
y1=[ 0  10];

y1=[y1 -fliplr(y1)];
x1=[x1 -fliplr(x1)];

quellfunc=@(s)interp1(x1,y1,s,'linear');
%%%%%%%%%%%%%%%%%%%%%%

[x,y]=meshgrid(linspace(-2,2,12),linspace(-2,2,12));

[m,n]=size(x);

phi=x.*geschw;
psi=y.*geschw;

for a=1:m
    for b=1:n
        if ((x(a,b)<=d && x(a,b)>=-d) && y(a,b)==0)
        psi(a,b)=nan;
        else
         phi(a,b)=phi(a,b)+1/(2*pi)*quad(@(s)quellfunc(s).*log(sqrt((x(a,b)-s).^2+y(a,b).^2)),-d,d);
         psi(a,b)=psi(a,b)+1/(2*pi)*quad(@(s)quellfunc(s).*atan2(y(a,b),(x(a,b)-s)),-d,d);
        end
    end
end

meshc(x,y,psi);