%fd example
%solve du/dt = d2u/dx2
%BCs du/dx = u  at x=0   du/dx = -u  at x=1
%Initial conditions u=1

clc;
clear all;

%parameters
dx = 0.1;
dt = 0.0025;
Nsteps = 40;

%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) = 1.0;
end

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

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

 %BCS (compute before updade u)
 u_init = u(1) + 2*r*(u(2) - (1+dx)*u(1));
 u_end = u(Nx) + 2*r*(u(Nx-1) - (1+dx)*u(Nx));
    
 %update u
 u = u + r*(circshift(u,[0 +1]) - 2*u + circshift(u,[0 -1]));
 u(1) = u_init;
 u(Nx) = u_end;

 %plot result
 plot(x,u);
end
