title: Applying neumann boundary conditions to the diffusion equation
body: <p>I've plotted a code for the the numerical solution to the diffusion equation du/dt=D(d^2 u/dx^2) + Cu where u is a function of x and t - I've solved it numerically and plotted it with the direchtlet boundary conditions u(-L/2,t)=u(L/2,t)=0, with the critical length being the value before the function blows up exponentially, which I have worked out to be pi. I'm trying to change the right boundary condition to be the Neumann boundary condition u_x(L/2)=0 which should decrease the critical length and cause the function to increase exponentially, but I'm not too sure exactly how to do this and mine doesn't seem to work quite right - could someone see if they could identify where I've gone wrong? Thanks!</p>
  
  <pre><code>L=np.pi # value chosen for the critical length
  s=101 # number of steps in x
  t=10002 # number of timesteps
  ds=L/(s-1) # step in x
  dt=0.0001 # time step
  D=1 # diffusion constant, set equal to 1
  C=1 # creation rate of neutrons, set equal to 1
  Alpha=(D*dt)/(ds*ds) # constant for diffusion term
  Beta=C*dt # constant for u term
  
  x = np.linspace(-L/2, 0, num=51)
  x = np.concatenate([x, np.linspace(x[-1] - x[-2], L/2, num=50)]) # setting x in the specified interval
  
  u=np.zeros(shape=(s,t))
  u[50,0]=1/ds # delta function
  for k in range(0,t-1):
      u[0,k]=0 #direchtlet boundary condition
      for i in range(1,s-1):
          u[i,k+1]=(1+Beta-2*Alpha)*u[i,k]+Alpha*u[i+1,k]+Alpha*u[i-1,k] # numerical solution 
      u[s-1,k+1]=u[s-2,k+1] # neumann boundary condition
      if k == 50 or k == 100 or k == 250 or k == 500 or k == 1000 or k == 10000: # plotting at times
          plt.plot(x,u[:,k])
  
  plt.show()
  </code></pre>