%% Heat equation u_t=u_xx - explicit finite difference scheme

function heat_explicit

% This is a fancy way to define a function (primary) which calls another 
% function (secondary), defined at the bottom of this file

% Copyleft 2010 Radu C. Cascaval, UCCS

%% time step and spatial step size
dt=0.0012;
dx=0.05;

% number of time iterations
N=50;

% final time of the computation
Tf=N*dt;

% the mesh ratio
mu=dt/(dx)^2;

tvals = 0:dt:Tf;
xvals = 0:dx:1;

u0 = @(x) 1-2*(x-1/2).*sign(x-1/2);
init = u0(xvals);

uvals = heatD(xvals, tvals, init, [0 0]);

surf(xvals, tvals, uvals') % plots the graph of u=u(x,t)
xlabel x; ylabel t; zlabel u  
title(strcat('mesh size dx=',num2str(dx),...
    ', ratio mu=dt/(dx)^2= ',num2str(mu)))



%% This is where the primary function ends and a secondary function starts

function u = heatD(x, t, init, bdry)
% solve the 1D heat equation on the rectangle described by
% vectors x and t with u(t(1), x) = init and Dirichlet boundary
% conditions u(t, x(1)) = bdry(1), u(t, x(end)) = bdry(2).

N=length(t);
J=length(x);

dt=mean(diff(t));
dx=mean(diff(x));

mu = dt/dx^2;

u = zeros(J,N);
u(:,1) = init;

for n = 1:N-1
	u(2:J-1,n+1) = mu*u(3:J,n) + (1-2*mu)*u(2:J-1,n)+mu*u(1:J-2,n);
	u(1,n+1) = bdry(1);
	u(J,n+1) = bdry(2);
end
