function diffusion_1d_SL
%% Solving the time dependent 1d diffusion problems using
%% separation of variables and a sine series, effectively a
%% Sturm-Liouville (SL) approach.
%% Greenberg 18.3.1, example 1: from eqn 1a,b,c to 28, pp 954-958.
%% Eli, APM105b

clf; clear all

%% choose initial conditions scenario:
for icase=0:5;  %% 0 to 5

  %% diffusion coeff (cm^2/sec):
  alpha=sqrt(1.14);
  %% size of domain (cm):
  L=10.0;
  %% side b.c.:
  if icase==0 || icase==1
    u1=0; u2=0;
  elseif icase==2 || icase==3
    u1=0; u2=75;
  elseif icase ==4 || icase ==5
    u1=50; u2=50;
  end
  %% number of Fourier coeffs:
  Nf=30; K=zeros(Nf,1);
  if Nf>30; 
    disp('quad doesn''t work well with large Nf, try Nf=<30.');
  end
  %% array with x-axis for plotting:
  x_axis=0:0.02:L;
  
  %% calculate fourier coeffs:
  TOL=1.e-9; %% tolerance for integral calculating fourier coeffs
  for n=1:Nf
    K(n)=(2/L)*quad(@(x) (F_ic(x,n,L,icase)-(u1+(u2-u1)*x/L)) .*sin(n*pi*x/L),0,L,TOL);
  end


  %% plot an animation of the solution:
  times=[0:0.001:0.05,0.1:0.1:20];
  for t=times
    %% Calculate the actual solution (22) of the 1d diffusion equation.
    uu=u1+(u2-u1)*x_axis/L;
    for n=1:Nf
      uu=uu+K(n)*sin(n*pi*x_axis/L)*exp(-(n*pi*alpha/L)^2*t);
    end
    plot(x_axis,uu);
    h1=xlabel('x');
    h2=ylabel('Temperature');
    h3=title(sprintf('t=%6.3f; solution of 1d diffusion eqn',t));
    set([gca,h1,h2,h3],'fontsize',14);
    axis([0 L 0 130]);
    %% plot actual initial conditions
    hold on
    plot(x_axis,F_ic(x_axis,n,L,icase),'r')
    hold off
    if t==0; pause; end
    pause(0.1)
  end

  pause
end %% end of icase loop

function out=F_ic(x,n,L,icase)
%% Specified initial conditions, needs to be able to accept a
%% vector x and return a vector with the i.c. at the corresponding
%% locations
if icase==0
  out=70*(exp(-(x-0.5*L).^2/(L*0.2)^2));
elseif icase==1 || icase==2
  out=x*0+100;
elseif icase==3
  out=70*(exp(-(x-0.35*L).^2/(L*0.2)^2)+exp(-(x-1*L).^2/(L*0.2)^2));
elseif icase==4
  out=50+20*sin(3*pi*x/L);
elseif icase==5
  out=50+20*sin(10*pi*x/L);
end
