%  Heat Diffusion in a Thin Insulated Wire

clear all; close all;
L = 1.0;                           %  Lenth of the Wire
tfinal = 1;                        %  Final Time

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

rho = 1.; cond = 1/(pi*pi); spheat = 1.0;    % Material parameters

dt = dx*dx*rho*spheat/(cond*2.1);
maxk = fix(tfinal/dt),              %  Number of Time Steps

b = dt/(dx*dx); a = cond/(spheat*rho); d = a*b;

%  Initial Temperature
for i = 1:n+1
        x(i) =(i-1)*dx;
        u(i,1) =sin(pi*x(i));
        error(i,1) =0;
end

for k=1:maxk+1                         %Boundary Temperature
        u(1,k) = 0;
        u(n+1,k) = 0;
        time(k) = (k-1)*dt;
	error(1,k) = 0;  error(n+1,k) = 0;
end

for k=1:maxk                % Time Loop
   for i=2:n;               % Space Loop 
      u(i,k+1) =0.*dt/(spheat*rho)+(1-2*d)*u(i,k) + d*(u(i-1,k)+u(i+1,k));
      error(i,k+1)= abs(u(i,k+1)-exp(-time(k))*sin(pi*x(i)));
   end
end

mesh(x,time,u');  figure(2); mesh(x,time,error')

format short e; max(error(:,maxk))