% Math 3150 Fall 2008
% Example 3.6.3 - heat equation with Robin BC

% first we get a good guess of the first few roots mu_n by looking at where
% the two following curves intersect
mu = linspace(0,5*pi,300);
plot(mu,tan(mu),mu,-mu);
axis([0 5*pi -15,5]);

amus = [3 5 8 11 14.3 17.5];
% then we let Matlab refine the roots by using fsolve
% (fzero tends to choke on this problem, so we call the optimization toolbox)
for i=1:length(amus),
 mus(i) = fsolve(@(mu) tan(mu) + mu,amus(i));
end;
mus

toprint = 0; % 1 create print quality plots, 0 animate
c=1;         

% space variable
x = linspace(0,1);

% time variable will take values in this array
ts=linspace(0,1,21); % 21 equally spaced points from 0 to 1, included
		     % depending on c you may need to change this


figure(2); clf;
syms y
if (toprint) hold on; end; % do not erase previous plots
for it=1:length(ts), % for every time
 t = ts(it);
 s = zeros(size(x)); %initialize partial sum to zero
 for n=1:length(mus),
  cn = int(y*(1-y)*sin(mus(n)*y),0,1)/int(sin(mus(n)*y)^2,0,1);
  cn = double(cn); % convert symbolic variables to doubles
  s = s + cn*exp(-mus(n)^2*t)*sin(mus(n)*x);
 end;
 plot(x,s,'b',x,x.*(1-x),'r'); % plot the partial sum and initial condition
 axis([0 1 0 0.25]); % adjust axis so that all plots are on the same scale
 pause; % wait for user input
end;
if (toprint) hold off; end; % erase previous plots
