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

% parameters for Newton's method
maxit=60; % maximum number of iterations
tol  =1e-6; % tolerance

% 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,6*pi,300);
figure(1); clf;
plot(mu,tan(mu),mu,-mu);
axis([0 6*pi -20,5]);


% these initial guesses work empirically
% and were determied from looking at the plot
% (the roots seem to be just after odd multiples of pi/2)
amus = (1:2:13)*pi/2 + 1e-3;
mus = zeros(size(amus));

% Use Newton's method to find the roots of tan(mu) + mu,
% i.e. the mu for which tan(mu) + mu = 0
% We use the initial guess amus(i)
for i=1:length(amus)
 % Newton's method
 m = amus(i);
 for k=1:maxit,
  % Newton's iteration is m = m - f(m)/f'(m)
  mnew = m - (tan(m) + m)/(2+tan(m)^2);
  % check convergence
  if (abs(mnew-m)/abs(mnew) < tol) 
   m=mnew; 
   break; 
  else
   m=mnew;
  end;
 end;
 mus(i) = m;
 fprintf('root #%2d close to %6.4g is %13.10g (in %2d iter)\n',i,amus(i),m,k);
end;

% add these roots to the plot just to make sure everything is right
hold on;
plot(mus,tan(mus),'r*');
hold off;

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

% space variable (used for integration as well)
N = 100;
h = 1/(N-1);
x = h*(0:N-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;
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),
  % compute dot product ( x*(1-x) , sin(mus(n)*x) ) 
  % (numerical integration using composite trapezoidal rule)
  fsin   = h*trapz(x.*(1-x).*sin(mus(n)*x));

  % compute dot product ( sin(mus(n)*x) , sin(mus(n)*x) )
  % (numerical integration using composite trapezoidal rule)
  sinsin = h*trapz(sin(mus(n)*x).*sin(mus(n)*x));

  % the coefficients are
  cn = fsin/sinsin;

  % compute sum
  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
