Help! Why is my code so slow/ does this even work properly? Using a crank-Nicholson scheme to solve a PDE.

Hello,

I am using a numerical scheme (crank-Nicholson scheme) to solve a PDE as part of an assignment.

This is an diffusion-advection-reaction problem.

The equation in question is the following.

dC/dt= D*[d^2*C/dx^2]- ad*[dC/dx]- R(x)

The first term and second terms have constant coefficients, while the reaction term is a function of x. Part of using the crank Nicholson scheme is inverting a matrix.As you can see, with a great number of points this can be a problem (when using the \ command). I tried using the conjugate gradient code provided by mathworks, but either 1) I am not using it correctly or 2) the way my code works makes everything slow. The code is below. Any tips or suggestions will be appreciated!

~MsCaptainFlan **************************************************************

function [pp ] = CrankNick3
%Solves the advection-diffusion-reaction PDE using the crank-nicholson
%scheme
% instead of using the command AA\BB, the conjugate graduent method is used
% to speed up computation
%the conjugate gradient code was obtained from MathWorks
xl=0;xr=1;%boundaries in space
yb=0;yt=1; %boundaries in time
M=50; %number of points in space
N=100; %number of points in time
dx=(xr-xl)/M;dt=(yt-yb)/N;            % step sizes
m=M-1; n=N;
x=xl:dx:xr;tt=yb:dt:yt;
xx=x(1:m);

G0=10; %uM
%coefficients 
global D W k b tol
D=8e-8; %diffusion
W=-1*.1e-2;%advection
tol=1e-5;

%Rxn should be a m*1 vector. The rxn coefficient of the PDE varies in space
k=.15; %Rcox from Burdige
b=.1; %the beta reaction term
Rxn=k*exp(-b*xx'); 
%Coefficients used to generate matrix AA and BB
global alpha beta gama
alpha=dt*D/(2*dx);
beta=-dt*W/(4*dx);
gama=-dt*Rxn./2;

global c1 c2 c3 c4 c5 c6
c1=-alpha-beta;
c2=1+2*alpha-gama(1:m);
c3=-alpha+beta;
c4=alpha+beta;
c5=1-2*alpha+gama(1:m); 
% c5=gama(1:m);
c6=alpha-beta;

%Generate AA matrix
% define tridiagonal matrix AA
AA=diag(c2.*ones(m,1))+diag(c3.*ones(m-1,1),1);
AA=AA+diag(c1*ones(m-1,1),-1); 
%the next few lines define the boundary conditions
AA(1,:)=0; %the first row of any colum must be zero
AA(1,1)=1; %the first row, first column value must be one
AA(m,:)=0; %the last row, any colum must be zero
AA(m,m)=1; %the last row, last column must be one

%form initial matrix B0 with some concentration profile
% B0=G0;
Bb=zeros(m,1); %this forms the initial b matix
Bb(1,1)=(G0); %this defines the initial conditions
w0=AA\Bb; 

BB=zeros(m,1);

BB=Bb; %for the first iteration of BB, let it be initial condition vector Bb
w=w0; %for the first iteration of w, let it be initial condition vector w0
tic
for j=1:n,
    for p=2:m,
        BB(p)=(c4*w(p))+(c5(p).*w(p))-(c6*w(p));

          BB;
          gcf; 
          xp=x';

           fignodiff=plot(xp(1:49,1),w);hold on;

          xlabel('depth');ylabel('c');
          drawnow;
      end
  w=conjgrad(AA,BB,tol);

end
toc
end

function w = conjgrad(AA,BB,tol)
% CONJGRAD  Conjugate Gradient Method.
%   X = CONJGRAD(A,B) attemps to solve the system of linear equations A*X=B
%   for X. The N-by-N coefficient matrix A must be symmetric and the right
%   hand side column vector B must have length N.
%
%   X = CONJGRAD(A,B,TOL) specifies the tolerance of the method. The
%   default is 1e-10.
%
% Example:
%{
  n = 6000;
  m = 8000;
  A = randn(n,m);
  A = A * A';
  b = randn(n,1);
  tic, x = conjgrad(A,b); toc
  norm(A*x-b)
%}
% By Yi Cao at Cranfield University, 18 December 2008.
%
    if nargin<3
        tol=1e-10;
    end
    r = BB + AA*BB;
    y = -r;
    z = AA*y;
    s = y'*z;
    t = (r'*y)/s;
    w = -BB + t*y;

      for k = 1:numel(BB);
         r = r - t*z;
         if( norm(r) < tol )
              return;
         end
         B = (r'*z)/s;
         y = -r + B*y;
         z = AA*y;
         s = y'*z;
         t = (r'*y)/s;
         w = w + t*y;
      end
   end