%This code plots solutions to diffusion
%in the infinite cylinder with radial flow
%variable temperature

%solve 
%du/dt = k LAPLACE(u)
%u(a,t) = 0 a = radius at boundary
%u(r,0) = f(r)

a = 1;
f = @(x) 10*(x-a).*(x-0.5);

% 5*sin(pi*x);

N = 25;
k = .001; %a^2
Nt = 20;
delt = 0.001;

r = [0:0.05:a];
J0 = @(x) besselj(0,x);
J1 = @(x) besselj(1,x);


%need N roots of J0
for n = 1:N
z(n) = fzero(J0,[(n-1)*pi, n*pi]);
end

% % x = 0:pi/50:N*pi;
% % y = J0(x);
% % plot(z,zeros(1,N),'o',x,y,'-')
mu = z/a;

%First make coefficents of f(r)
for n = 1:N
F = @(x) x.*f(x).*J0(mu(n).*x);
tmp= quad(F,0,a);
A(n) = 2*tmp/(J1(a*mu(n)))^2;
end


[th,r2] = meshgrid((0:5:360)*pi/180,r);
[X,Y] = pol2cart(th,r2);

%Compute u
u = zeros(Nt,length(r));
for i = 1:Nt
    for j = 1:length(r)
        for n = 1:N
            u(i,j) = u(i,j) + A(n)*exp(-k*mu(n)^2*delt*i)*J0(mu(n)*r(j));
        end
    end
    
    for k = 1:max(size(X))
       U(:,k,i) = u(i,:)';
    end
end
for j=1:Nt
  surf(X,Y,U(:,:,j))
  axis([-1 1 -1 1 0 5] )
  pause
end
  