% parameterize space
dx = 1;
dt = .1;
x = -10:dx:10;
t = 0:dt:10;

% starting values C(x,t=0) = n0*delta(x)
% dc(x,t) = D*d^2c(x,t)/dx^2
n0 = 5;
D = 2;

% set size of c(x,t) to be filled in later
c = NaN*ones(length(x),length(t));

% set initial conditions: C(x,t=0) = n0*delta(x)
c(:,1) = zeros(length(x),1);
x0_ind = ceil(length(x)/2);
c(x0_ind,1) = n0/dx;

% run equation:
for it = 1:length(t)-1 %at each time step do the following:
    dcdx = diff(c(:,it))/dx; % compute dc/dx at each time step
    dcdx2 = [0;diff(dcdx);0]/dx; % compute d(dc/dx)/dx at each time step
    dcdt = D*dcdx2; % Diffustion equation
    c(:,it+1) = c(:,it) + dcdt*dt; % update c
end

% plot
figure(1)
[T,X] = meshgrid(t,x);
surf(T,X,c);
xlabel('t');
ylabel('x');
zlabel('c');

% Analytic solution
c_analytic(:,1) = c(:,1);
for ix = 1:length(x);
    for it = 2:length(t);
        c_analytic(ix,it) = n0/sqrt(4*pi*D*t(it))*exp(-(x(ix))^2/(4*D*t(it)));
    end
end
figure(2)
surf(T,X,c_analytic);
xlabel('t');
ylabel('x');
zlabel('c');