function [xsource,msource]=AxiSymmBody(L,c,epsilon,N)
%Give the location of the sources.
xsource=zeros(N,1);
for j=1:N
    xsource(j)=-L+2*L*j/(N+1);
end
%Give the body shape  fbody ,  and its derivative fprime.
fbody=zeros(N,1);
fprime=zeros(N,1);
for j=1:N
    fbody(j)=bodyshape(xsource(j),L,epsilon);
    fprime(j)=bodyprime(xsource(j),L,epsilon);
end
%Set up the coefficient matrix  Acoef .
Acoef=zeros(N,N);
for i=1:N
    for j=1:N
        diffij=xsource(i)-xsource(j);
        Acoef(i,j)=(-diffij*fprime(i)+fbody(i))/(diffij^2+(fbody(i))^2)^(3/2);
    end
end
%Set up the right-hand side vector  brhs .
brhs=zeros(N,1);
for i=1:N
    brhs(i)=4*pi*c*fprime(i);
end
%Solve the matrix equation.
msource=Acoef\brhs;
figure(1)
plot(xsource,msource)
%plot the body and the velocity vector.
oneeL=1.375*L;
stepL=0.25*L;
[X,Z]=meshgrid(-oneeL:stepL:oneeL,-oneeL:stepL:oneeL);
figure(2)
plot(xsource,fbody,'r')
hold on
fbodymin=-fbody;
plot(xsource,fbodymin,'r')
U=zeros(size(X));
W=zeros(size(X));
for j=1:N
    U=U+msource(j)*(X-xsource(j))./((X-xsource(j)).^2+Z.^2).^(3/2);
    W=W+msource(j)*Z./((X-xsource(j)).^2+Z.^2).^(3/2);
end
U=U/(4*pi);
W=W/(4*pi);
U=U+c;
quiver(X,Z,U,W)
hold off
function f=bodyshape(x,L,epsilon)
f=epsilon*(L^2-x^2)/(L^2);
function fdash=bodyprime(x,L,epsilon)
fdash=-2*epsilon*x/(L^2);