
% Steady state temperature of a semi-infinite cylinder of radius=1
% and bottom temp=100 and side temp=0

format compact;
Temp = 100;			% temperature of bottom of cylinder
zfin = 1.5;			% plot up to z=zfin

NptsR=101;
NptsZ=51;
Nplots=25;
Nrows=5;			% number of columns
Ncols=Nplots/Nrows;
Nskip=1;			% use every Nskip coefs
r=linspace(0,1,NptsR);
z=linspace(0,zfin,NptsZ);
[R,Z]=meshgrid(r,z);

% the first 201 zeros of J_0 (computed with Maple's evalf(BesselJZeros(0;1..201));)
BesselJ0zeros=[2.404825558; 		
5.520078110; 8.653727913; 11.79153444; 14.93091771; 18.07106397; 21.21163663; 24.35247153; 27.49347913;
30.63460647; 33.77582021; 36.91709835; 40.05842576; 43.19979171; 46.34118837; 49.48260990; 52.62405184;
55.76551076; 58.90698393; 62.04846919; 65.18996480; 68.33146933; 71.47298160; 74.61450064; 77.75602563;
80.89755587; 84.03909078; 87.18062984; 90.32217264; 93.46371878; 96.60526795; 99.74681986; 102.8883743;
106.0299309; 109.1714897; 112.3130503; 115.4546127; 118.5961766; 121.7377421; 124.8793089; 128.0208770;
131.1624463; 134.3040166; 137.4455880; 140.5871604; 143.7287336; 146.8703076; 150.0118825; 153.1534580;
156.2950343; 159.4366112; 162.5781887; 165.7197667; 168.8613454; 172.0029245; 175.1445041; 178.2860842;
181.4276647; 184.5692456; 187.7108270; 190.8524087; 193.9939907; 197.1355731; 200.2771558; 203.4187388;
206.5603221; 209.7019057; 212.8434896; 215.9850737; 219.1266580; 222.2682426; 225.4098274; 228.5514125;
231.6929977; 234.8345831; 237.9761688; 241.1177546; 244.2593406; 247.4009267; 250.5425130; 253.6840995;
256.8256861; 259.9672729; 263.1088598; 266.2504469; 269.3920341; 272.5336214; 275.6752088; 278.8167963;
281.9583840; 285.0999718; 288.2415596; 291.3831476; 294.5247357; 297.6663239; 300.8079121; 303.9495005;
307.0910889; 310.2326775; 313.3742661; 316.5158548; 319.6574435; 322.7990324; 325.9406213; 329.0822103;
332.2237994; 335.3653885; 338.5069777; 341.6485669; 344.7901563; 347.9317456; 351.0733351; 354.2149246;
357.3565141; 360.4981037; 363.6396934; 366.7812831; 369.9228729; 373.0644627; 376.2060525; 379.3476424;
382.4892324; 385.6308224; 388.7724124; 391.9140025; 395.0555926; 398.1971828; 401.3387730; 404.4803632;
407.6219535; 410.7635438; 413.9051341; 417.0467245; 420.1883149; 423.3299053; 426.4714958; 429.6130863;
432.7546769; 435.8962675; 439.0378581; 442.1794487; 445.3210393; 448.4626300; 451.6042207; 454.7458115;
457.8874023; 461.0289930; 464.1705839; 467.3121747; 470.4537656; 473.5953565; 476.7369474; 479.8785383;
483.0201293; 486.1617203; 489.3033113; 492.4449023; 495.5864933; 498.7280844; 501.8696755; 505.0112666;
508.1528577; 511.2944488; 514.4360400; 517.5776312; 520.7192224; 523.8608136; 527.0024048; 530.1439961;
533.2855873; 536.4271786; 539.5687699; 542.7103612; 545.8519526; 548.9935439; 552.1351353; 555.2767266;
558.4183180; 561.5599094; 564.7015008; 567.8430923; 570.9846837; 574.1262752; 577.2678666; 580.4094581;
583.5510496; 586.6926411; 589.8342326; 592.9758242; 596.1174157; 599.2590073; 602.4005988; 605.5421904;
608.6837820; 611.8253736; 614.9669652; 618.1085568; 621.2501485; 624.3917401; 627.5333317; 630.6749234];


b=zeros(1,Nplots);

% Fourier coefs for U0(x) temperature
for jj=1:Nskip:Nplots
  b(jj) = (2*Temp)/(BesselJ0zeros(jj)*besselj(1,BesselJ0zeros(jj)));
end

b(1:Nskip:Nplots)

U=zeros(NptsZ,NptsR);			% initial steady state distribution (Fourier)

fs=8;
figure(3)
clf
plotsols=1;
for jj=1:Nskip:Nplots		% build Fourier series
  Un=b(jj)*exp(-BesselJ0zeros(jj)*Z).*besselj(0,BesselJ0zeros(jj)*R);
  U=U+Un;
  if(plotsols==1)
   subplot(Nrows,Ncols,jj)
   set(gca,'FontSize',[fs]);
   surf(R,Z,Un);
   shading interp
   zl = ['u(b_{',num2str(jj),'})'];
   xlabel('r','Fontsize',fs); ylabel('z','Fontsize',fs); 
   view([30 60])
   set(gca,'FontSize',[fs]);
   if(jj==1)
     axis tight
     ax=axis
     ax(5)=-ax(6); 
   end
   axis(ax);
   axis tight;
   ax=axis;
   text((ax(2)-ax(1))/2,ax(4),ax(6),zl,'Fontsize',fs);
   pause(0.0001)
  end
end

% Now plot the Fourier combo:
fs=16;
figure(4)
clf

  subplot(2,1,1)
  set(gca,'FontSize',[fs]);
  surf(R,Z,U);
  shading interp
  xlabel('r','Fontsize',fs); ylabel('z','Fontsize',fs); zlabel('z','Fontsize',fs);axis tight
  view([60 32])
  set(gca,'FontSize',[fs]);

  subplot(2,1,2)
  set(gca,'FontSize',[fs]);
  surf(R,Z,U);
  shading interp
  colormap(jet)
  view([0 90])
  xlabel('r','Fontsize',fs); ylabel('z','Fontsize',fs); zlabel('z','Fontsize',fs);axis tight
  colorbar
  set(gca,'FontSize',[fs]);




