%fd example
%solve du/dt = d2u/dx2
%BCs u(0) = u(1) = 0
%Initial conditions u=2x (x<0.5) & u=2(1-x) (x>0.5)

clc;
clear all;

%parameters
dx = 0.1;
dt = 0.001;
Nsteps = 20;

%resulting factor
r = dt/(dx^2);

%number of solution points
Nx = floor(1.0/dx) + 1;

%vector of solution point locations
x = (0.0:1.0:Nx-1)*dx;

%initialize solution (at first time step)
for i=1:Nx
 u(i) = 2*x(i)*( x(i) <= 0.5 ) + 2*(1-x(i))*( x(i) > 0.5 );
end

%plot initial values
plot(x,u);
xlabel('x');
ylabel('u');
title('Smith Ex. 2.1');
hold on;

%time loop
for itime = 1:Nsteps
 itime*dt

 %update u
 u = u + r*(circshift(u,[0 +1]) - 2*u + circshift(u,[0 -1]));

 %BCS
 u(1) = 0.0;
 u(Nx) = 0.0;

 %plot result
 plot(x,u);
end
