function [ u ] = reaction_diffusion( R, D, u0, num_t, dt)
%
% reaction_diffusion.m, (c) M Roughan, 2018
%
% created: 	Wed Apr 11 2018 
% author:  	M Roughan 
% email:   	matthew.roughan@adelaide.edu.au
% 
% reaction-diffusion system on a 1D ring
%    du/dt = R(u) + D \nabla^2 u
%
% INPUTS:
%       R      = function governing the reaction rates
%       D      = diffusion coefficients (should be a MxM diagonal array)         
%       u0     = starting concentrations ( MxN array, M=no of reagents, N = no. of cells)
%       num_t  = number of time steps
%       dt     = size of time steps
%         
% OUTPUTS:        
%       u      = reagent concentrations after num_t steps  
%  
% 
D = diag(diag(D)); % make D diagonal if it isn't

num_x = size(u0,2); % number of cells
M = size(u0,1); % number of reagents
dx = 1/num_x;
u = u0;
u_new = zeros(size(u0)); % initialise everything to zero

for j=1:num_t
  for i=2:num_x-1  
    u_new(:,i) = u(:,i) + dt*R(u(:,i)) + ...
	D*dt*( u(:,i+1) - 2*u(:,i) + u(:,i-1) )/dx^2; 
  end 
  % circular edges ...
  u_new(:,1) = u(:,1) + dt*R(u(:,1)) + ...
      D*dt*( u(:,2) - 2*u(:,1) + u(:,end) )/dx^2; 
  u_new(:,end) = u(:,end) + dt*R(u(:,end)) + ...
      D*dt*( u(:,1) - 2*u(:,end) + u(:,end-1) )/dx^2; 
  
  if any(any(u_new < 0))
    % if we detect instability then stop
    j
    break
  end
  u = u_new;
end
