%  Heat Diffusion in a Thin Insulated Wire

clear all; close all;
global beta c

beta = 2;   c = 0.5;               %  material parameters
L = 5.0;                           %  Length of the Wire
tfinal = 1;                        %  Final Time

n = 40.; h = L/n;                 %Number of Space Steps

dt = h;
maxk = fix(tfinal/dt),              %  Number of Time Steps

%  Initial Temperature
for i = 1:n+1
        x(i) = (i-1)*h;
        u0(i) = uexact(x(i),0);
end

plot(x,u0); hold

t = 0;

A = sparse(n+1,n+1); b=zeros(n+1,1); 
A(1,1) = 1; 
for i=2:n
  A(i,i) = 1/dt + beta/(h*h) + c/2;
  A(i,i-1) = -beta/(2*h*h);       A(i,i+1) = -beta/(2*h*h);
end
A(n+1,n+1) = 1;


for k=1:maxk                % Time Loop

   b(1) = 0;   b(n+1) = 0;      % The boundary condition
   for i=2:n;               % Space Loop 
      b(i) = f(x(i),t+0.5*dt) + u0(i)/dt+ 0.5*beta*( u0(i-1) -2*u0(i) + ...
             u0(i+1) )/(h*h) -c*u0(i)/2;
   end
   u1 = A\b;          % Solve the system of equations.
   
   u0 = u1;   plot(x,u1)
   t = t + dt;
end

for i=1:n+1
  err(i) = u1(i)- uexact(x(i),t);
end

  figure(2); plot(x,err)
  format short e
  e2 = [norm(err,1)/n,norm(err,2)/sqrt(n),norm(err,inf)],