% Script for solving heat equation of HW1, prob. 1g ,
%   dpsi/dt = a* d^2psi/dx^2
% with IC              psi(x,0) = 1 - cos(2*pi*t),
%      Neumann BCs dpsi/dx(0,t) = dpsi/dx(1,t) = 0.
% We use FTCS method with grid spacing dx = 1/N and image points at 
% x(1) = -1/N, (N+1)/N:
%   d_t^F phi_j^n = a*(d_x)^2 phi_j^n
% with 
%   x_j = (j-2)/N, j = 1,...,N+3
% and BCs phi_1^(n+1) = phi_3^(n+1) and phi_(N+3)^(n+1) = phi_(N+1)^(n+1)
% The script plots the solution and error for N = 5, 10, 20 and 40.

  a = 1;
  nu = 0.4;
  T = 0.016;
  L = 1;
  
  mrange = 1:4;   % Vector of indices for looping over resolution
  label = {'b+','gx','r--','k-'}; % Plot labels
  Nrange = zeros(4); % Initialize sequentially filled arrays
  maxerr = zeros(4);

  for m = mrange;
    N = 5*2^(m-1);
    Nrange(m) = N; % Used later for plotting
    dx = L/N; 
    x = dx*(-1:(N+1)); % Range of x includes image points
    phi = 1 - cos(2*pi*x); % IC (note this automatically satisfies BCs)
    dt = nu*dx^2/a;
    nt = round(T/dt);
    for n = 1:nt
      phi(2:(N+2)) = phi(2:(N+2))...
	  + nu*(phi(3:(N+3))-2*phi(2:(N+2))+phi(1:(N+1))); % Update interior
      phi(1) = phi(3);  % Update Neumann BCs
      phi(N+3) = phi(1);
    end
    if(m==1)

      % Set up plot of solution (except image points) for the various N's.

      subplot(2,2,1)
      plot(x(2:(N+2)),phi(2:(N+2)),label{1})
      xlabel('x')
      ylabel('phi(x,T)')
      hold on
    else
      plot(x(2:(N+2)),phi(2:(N+2)),label{m})
    end
    phiex = 1 - cos(2*pi*x)*exp(-4*pi^2*T);
    maxerr(m) = max(phi - phiex);
  end
  legend('N=5','N=10','N=20','N=40')
  title(['FTCS solutions: \nu = 0.4, T = ' num2str(T)])
  hold off
  
  subplot(2,2,2)
  loglog(Nrange,maxerr,'x',Nrange,maxerr(4)*(Nrange/Nrange(4)).^(-2),'--')
  xlabel('N')
  ylabel('Max-norm error')
  legend('FTCS','O(N^{-2})')
  title('FTCS error convergence (dx = N^{-1})')