function [x,t,u,Frames] = heat1d
% Heat diffusion in a wire

% Setup
L = 1.0;
T = 150;
K = 0.001;
n = 150;
dt = T/n;
l =20;
dx = L/l;

u = zeros(l+1,n+1);

% Create mesh
t = (0:n)*dt;
x = (0:l)*dx;

% Set i.c. and b.c.s
u(:,1) = u0(x)';

u(1,:) = bcLeft;

u(l+1,:) = bcRight;

alpha = K*dt/dx^2;

for k=1:n
    for i=2:l
        u(i,k+1) = alpha*u(i-1,k) + (1 - 2*alpha)*u(i,k) + alpha*u(i+1,k)+ dt*f;
    end
end

for k=1:n+1
    plot(x,u(:,k));
    axis([0 1 0 10]);
    Frames(k) = getframe;
end

% Some functions

function f = f
f = 0.1;
%f = 0.0;

function u0 = u0(x)
u0 = 0.0;
% for i = 1:length(x)
%    if x(i) < 0.5
%        u0(i) = 5.0;
%   else
%        u0(i) = 0.0;
%    end
% end

function bcLeft = bcLeft
bcLeft = 0.0;

function bcRight = bcRight
bcRight = 0.0;
