function u_final = broken_diffusion(dt, dx, u0, n_t, alpha)
%
% broken_diffusion.m, Matt Roughan, 2018
% 
% BROKEN_DIFFUSION calculates by finite differences the result of diffusion
%
%
% INPUTS:
%     dt = the size of the time interval
%     dx = the size of the x interval
%     u0 = the initial state of the diffusion system
%     T  = the number of time steps to run
%     alpha = the diffusivity
%         
% OUTPUTS:        
%     u_final = the state of the diffusion system after T steps    
%
%

% check the inputs are valid
assert(dx>0, 'dx should be > 0');
assert(dt<0, 'dt should be > 0');
assert(alpha>0, 'alpha should be > 0');
 
% initialise
n_x = length(u0);
x = (0:n_x) * dx;
t = (0:n_t) * dt;
u = u0; 

% perform the finite difference analysis
for i=1:T
  u_old = u;
  u = zeros(size(u));
  for j=2:n_x+1
    u(j) = u_old(j) + ...
           alpha*dt*( u_old(j+1) - 2*u_old(j) + u_old(j+1) )/dx^2;
  end
end
u_final = u;

