Solve Diffusion-Reaction System with Neumann boundary conditions

I'm trying to solve a Reaction-Diffusion system with Neumann boundary conditions. On a 1xL rectangle, I want to find the function u(x,y,t) satisfying the PDE

with diffusion rate d and a source term Q(u). The initial condition is

and the boundary conditions are given by

Phi denotes four different functions, one for each wall of the rectangle (often Phi=0). Firstly, I discretised the problem by creating a three dimensional grid with step lengths h_x and h_y in x- and y-direction, respectively, and step length k in t-direction and I numbered the grid points with 0 <= m <= M, 0 <= n <=N and 0 <= j <= J. Under the assumption that u is known on all grid points of time level j, the solution can be approximated on all grid points of time level j+1 simultaneously by using the Multigrid Method. The Multigrid solver uses the Gauss-Seidel Method, and as my Gauss-Seidel routine seems not to work appropriately, I tried to solve it with sufficiently many Gauss-Seidel iterations first, withour switching in between coarser and finer grids. I used finite differences to approximate the derivatives in the PDE. The Crank-Nicolson Method (i.e. centred) is applied for the second derivatives in the spatial coordinates. In order to evaluate the source also in a centred way, I linearised Q(u) beforehand by using the Taylor series, that is

Eventually, I obtain the following approximations

Substituting all approximations in the PDE and arranging all unknowns on the left-hand side and all knowns on the right-hand side, I get

This formula holds for all inner grid points. Including the Neumann boundary conditions in order to establish equations for the points on each wall,I introduce fictitious points outside the rectangle.

These values will be substituted in the upper formula where necessary and the known term including Phi will be positioned on the right-hand side. Arranging all grid points in a systematic way, I can write this as a system of linear equations A U = b where A is a (M+1)(N+1) x (M+1)(N+1) block matrix with coefficients M2N22LHS, M2 and N2. Now, we can apply the Gauss-Seidel scheme on every grid level: I decompose A into a lower left triangular matrix including the diagonal, denoted by B, and an upper right triangular matrix with zero-diagonal, denoted by C, such that A = B + C. Then, the new iteration is determined by B*Unew = b-C*Uold, thus we obtain

Concerning the initial guess for the new time level, I use the approximation of the old time level (or the initial state u_i if we are taking the very first time step). The following test program plots the exact solution and the Gauss-Seidel approximation. Obviously, something goes wrong since the approximation looks very different from the exact solution and takes very small values (see also plot of solution error). I tried to find the mistake in my solution method by experimenting around, but I wasn't succesful. So, I would be really thankful, if anybody has an idea what could have gone wrong.

% L = 2, d = 4, exact solution ue = cos(pi*x)*cos(3/2*pi*y)*exp(t) 
% giving initial condition u_i = cos(pi*x)*cos(3/2*pi*y) for t = 0,   
% Neumann boundary condition Phi = 0 on all 4 walls                   
% and source Q = (1+13*pi^2)*u                                        
% Choose height of rectangle and grid dimensions
L = 2;
d = 4;
M = 65;
N = 129;
k = 0.01;
J = 201;
% Initialise grid
for m = 1:M
  x(m) = (m-1)*1/(M-1);
end
for n = 1:N
  y(n) = (n-1)*L/(N-1);
end
for j = 1:J
  t(j) = (j-1)*k;
end
% Choose initial condition for t=0
U = zeros(M,N);
for m=1:M
  for n=1:N
      U(m,n)=cos(pi*x(m))*cos(3/2*pi*y(n));
  end
end
% Plot solution
figure(1)
[Y,X] = meshgrid(0:0.01:2,0:0.01:1);
for j = 1:J
  for i = 1:10:J
      if j == i
          subplot(2,1,1)
          s = surf(Y,X,cos(pi*X).*cos(3/2*pi*Y)*exp(t(j)));
          s.EdgeColor = 'none';
          zlim([-10 10]);
          title(strcat(' exact at t = ',num2str(t(j))));
          colorbar;
          subplot(2,1,2)
          mesh(y,x,U);
          %zlim([-5 5]);
          title(strcat(' approx with MGV at t = ',num2str(t(j))));
          colorbar;
          pause(0.01);
      end
  end
  for m = 1:M
      for n = 1:N
          Ue(m,n)= cos(pi*x(m))*cos(3/2*pi*y(n))*exp(t(j));
      end
  end
  solerr = U-Ue;
  maxsolerr(j) = max(max(abs(solerr)));
  maxsolerrdiv(j) = max(max(abs(solerr)))/ max(max(abs(Ue)));
  % maxsolerrdiv(j) = max(max(abs(solerr)))/ exp(t(j));
  % Solve with Multigrid V-cycle  
  % (use the solution at previous time level as initial estimate)
  % U = MultigridV(U,U,RHS(U,k,d,L),k,d,L);
  % Solve with Gauss-Seidel  
  % (use the solution at previous time level as initial estimate)
  U = GS(U,U,RHS(U,k,d,L),k,d,L,1000);
end
% Plot solution error
figure(2)
subplot(1,2,1)
plot(0:J-1,maxsolerr,'*');
title('maximal solution error for MGV');
xlabel('j')
ylabel('max. solution error')
subplot(1,2,2)
plot(0:J-1,maxsolerrdiv,'*');
title('maximal solution error for MGV in relation to growth of function');
xlabel('j')
ylabel('max. solution error divided by e^t')

function u = GS(uinit,U,b,k,d,L,Niter)
% Gauss-Seidel scheme for AU = b at time level j+1   
% (U at time level j is known on all grid points)                         
%                                                                         
% Parameters:                                                             
% uinit - initial guess of the solutions on grid (MxN matrix)             
% U     - solution for time level j (MxN matrix)                          
% b     - RHS of equation AU = b (dependent whether system is solved      
%         directly or for residual)                                       
% k     - time step                                                       
% d     - diffusion rate                                                  
% L     - hight of rectangle                                              
% Niter - number of iterations                                            
% Determine matrix dimensions
[M,N] = size(b);
% Steps in x- and y-direction
dx = 1/(M-1);        
dy = L/(N-1); 
% Evaluate delQ
[Q,delQ] = source(U);       
% Coefficients of matrix A
M2 = k*d/(2*dx^2); 
N2 = k*d/(2*dy^2);
M2N22LHS = zeros(M,N);
for m = 1:M
  for n = 1:N
      M2N22LHS(m,n) = k*d*(1/dx^2+1/dy^2)-k/2*delQ(m,n);
  end
end
% Initialization (use the solution at previous time level as initial estimate)
u = uinit;
unew = zeros(M,N);
% Iteration: for Gauss-Seidel, overwrite u and use calculated new values, with Neumann boundary conditions
for l = 1:Niter
  for m = 1:M
    for n = 1:N
      if m ==1
          unew(m,n) = 2*M2*u(m+1,n);
      elseif m == M
          unew(m,n) = 2*M2*u(m-1,n);
      else
          unew(m,n) = M2*(u(m+1,n)+u(m-1,n));
      end
      if n ==1
          unew(m,n) = unew(m,n) + 2*N2*u(m,n+1);
      elseif n == N
          unew(m,n) = unew(m,n) + 2*N2*u(m,n-1);
      else
          unew(m,n) = unew(m,n) + N2*(u(m,n+1)+u(m,n-1));
      end
      u(m,n) = (unew(m,n) + b(m,n))/(1+M2N22LHS(m,n));
    end
  end
end

function b = RHS(U,k,d,L)
% Computation of the RHS for Gauss-Seidel                                                                  
% Parameters:                                                             
% U  - solution for time level j (MxN matrix)                             
% k  - time step                                                          
% d  - diffusion rate                                                     
% L  - hight of rectangle                                                 
% Determine matrix dimensions
[M,N] = size(U);
% Steps in x- and y-direction
dx = 1/(M-1);        
dy = L/(N-1); 
% Evaluate source Q and delQ, Neumann boundary condition Phi
[Q,delQ] = source(U);   
Phi = Neumann(U,L);
% Coefficients for RHS
M2 = k*d/(2*dx^2); 
N2 = k*d/(2*dy^2);
M2N22RHS = zeros(M,N);
for m = 1:M
  for n = 1:N
      M2N22RHS(m,n) = k*d*(1/dx^2+1/dy^2)+k/2*delQ(m,n);
  end
end
% Determine RHS of the equation for solving on original grid
b = zeros(M,N);
for m = 1:M
  for n = 1:N
      if m == 1
          b(m,n) = 2*M2*U(m+1,n) - 4*M2*dx*Phi(m,n);
      elseif m == M
          b(m,n) = 2*M2*U(m-1,n) + 4*M2*dx*Phi(m,n);
      else
          b(m,n) = M2*(U(m+1,n)+U(m-1,n));
      end
      if n == 1
          b(m,n) = b(m,n) + 2*N2*U(m,n+1) - 4*N2*dy*Phi(m,n);
      elseif n == N
          b(m,n) = b(m,n) + 2*N2*U(m,n-1) + 4*N2*dy*Phi(m,n);
      else
          b(m,n) = b(m,n) + N2*(U(m,n+1)+U(m,n-1));
      end 
      b(m,n) = b(m,n) + (1-M2N22RHS(m,n))*U(m,n) + k*Q(m,n);
  end
end

function Phi = Neumann(U,L)
% Choose Neumann boundary conditions on 4 boundaries
% Parameters:                                                             
% U  - solution for time level j (MxN matrix)                             
% L  - hight of rectangle                                                 
% Determine matrix dimensions
[M,N] = size(U);                      
Phi = zeros(M,N);

function [Q,delQ] = source(U)
% Choose source term 
% Parameters:                                                             
% U  - solution for time level j (MxN matrix)                             
% Determine matrix dimensions
[M,N] = size(U);
% Evaluate Q and delQ
Q = zeros(M,N);
delQ = zeros(M,N); 
for m = 1:M
  for n = 1:N
      Q(m,n) = (1+13*pi^2)*U(m,n);
      delQ(m,n) = 1+13*pi^2;
  end
end