function [] = heat()
% heat.m
%
% MATLAB function which runs a TFSC heat equation solver.
%
% DPN 2/15/12

clear all;
clf;

a = 0;
b = pi;
T = 1.0;
A = 0;
B = 0;

JJ = [10 20 40 80];

beta0 = 0.25;

fprintf('-----------\n');
fprintf('heat (tfsc)\n');
fprintf('-----------\n');
fprintf('a = %g  b = %g  T = %g\n',a,b,T);

for J=1:length(JJ)
  N(J) = floor(T*JJ(J)^2/(beta0*(b-a)^2));

  h(J) = (b-a)/JJ(J);
  k(J) = T/N(J);
  x = [a:h(J):b]';

  fprintf('J=%d N=%d h=%g k=%g\n',JJ(J),N(J),h(J),k(J));

  u0 = uex(x,0,a,b,1);
  uapprox = tfsc_heat(a,b,JJ(J),T,N(J),u0,A,B);
  uexact = uex(x,T,a,b,1);

  plot(x,uapprox,'b-o',x,uexact,'g-*');
  xlabel('x'); ylabel('u(x,T)');
  legend('approx','exact');
  pause(0.1);

  diff = uapprox-uexact;

  err(J) = norm(diff,inf)/norm(uexact,inf);
  fprintf('  rel error = %g\n',err(J));
end

c = polyfit(log(h),log(err),1);
fprintf('Rate of convergence: %g\n',c(1));

function [unew] = tfsc_heat(a,b,J,T,N,u0,A,B)
% tfsc_heat
%
% MATLAB function to use time-forward/space-centered finite difference
% scheme to simulate solutions of the heat equation:
%
% u_t = u_xx
% u(x,0) = u0(x)
% u(a,t) = A
% u(b,t) = B
%
% DPN 4/17/09

h = (b-a)/J;
k = T/N;
beta = k/h^2;

uold = u0;
unew = 0*uold;

for n=1:N
  t = n*k;
  unew(0+1) = A;
  for j=1:J-1
    unew(j+1) = beta*uold(j-1+1) + (1-2*beta)*uold(j+1)...
      + beta*uold(j+1+1);
  end
  unew(J+1) = B;

  uold = unew;
end

function [u] = uex(x,t,a,b,n)

omega = (2*pi*n/(b-a));

u = exp(-omega^2*t)*sin(omega*x);