% Solving the problem:
% u_t = u_xx  for 0 < x < 1 and 0 < t < 0.3
% u(x,0) = cos(Pi x)
% u(0,t) =   exp(-Pi^2 t)
% u(1,t) = - exp(-Pi^2 t)

function heat_02
clear all;
close all;

% setting the grid in space;
N = 10;
dt = 0.01;
Tmax = 0.3;
Nt = ceil(Tmax / dt);

% task (a)
% calculating and plotting the solution using the CN method;
U = Solve_Heat( N, dt, Nt, 0.5 );

x = linspace(0.0,1.0,N+1);
t = linspace(0.0,Tmax,Nt+1);
[T,X] = meshgrid(t,x);
figure(1); surf(X,T,U); 
xlabel('x'); ylabel('t'); zlabel('u'); 

% the analytical solution;
Uex = Sol_Exact( x, t );

figure(3); surf(X,T,Uex); 
xlabel('x'); ylabel('t'); zlabel('u'); 

% task (b)
% computing the errors corresponding to different values of theta;
dw = 0.125;
k = 1;
Nt0 = 5;
for w=0.0:dw:1.0  
  Uw = Solve_Heat( N, dt, Nt, w );
  Err(k,1) = w;
  Err(k,2) = max( abs( Uw(:,Nt0) - Uex(:,Nt0) ) );
  k = k + 1;
end;  
figure(2); plot(Err(:,1),Err(:,2),'-x');



return;


% The function solves the heat equation using the theta method;
function [u] = Solve_Heat( nx, dt, nt, theta )
  
global a;

  h = 1.0 / nx;
  x = 0 : h : 1;
  r = dt / h^2;
  
  u(2:nx,1:nt+1) = 0.0;
  u(1,1:nt+1)    = exp( - pi^2 * (0:nt) * dt );
  u(nx+1,1:nt+1) = - exp( - pi^2  * (0:nt) * dt );
  
% the initial condition and the RHS;
  u(:,1) = cos(pi*x');  
    
% assembling the matrices;
  A = (1 + 2.0 * r * theta)*diag(ones(1,nx-1)) - ...
      r * theta * (diag(ones(1,nx-2),1)+diag(ones(1,nx-2),-1));
  B = (1 - 2.0 * r * (1.0 - theta))*diag(ones(1,nx-1)) + ... 
      r * (1.0 - theta) * (diag(ones(1,nx-2),1)+diag(ones(1,nx-2),-1));  


  for k = 1 : nt
% the time-stepping solution
    G(1:nx-1,1) = B * u(2:nx,k);
    G(1)   = G(1) + r * ( (1.0 - theta) * u(1,k) + ...
			 theta * u(1,k+1) );
    G(nx-1) = G(nx-1) + r * ( (1.0 - theta) * u(nx+1,k) + ...
			 theta * u(nx+1,k+1) );
    u(2:nx,k+1) = A \ G;
    
  end;

  

% The function calculates the exact solution

function [u] = Sol_Exact( x, t )

  u = cos(pi * x)' * exp( - pi^2 * t );
 
  
  


