function [ML,MR] = matrixdefineTuring(N,lambda,x,dx,dim)
%%
%% This function sets up the matrices needed when using
%% the Crank-Nicolson method with derivative boundary
%% conditions at both x=0 and x=L.  
%%
%% (N+1)*dx = L
%%     xj = (j-1)*dx so j=1 corresponds to x=0
%%                  and j=N+2 corresponds to x=L.
%% There are N+2 unknowns (at j=1:N+2), so matrices
%%     need to be (N+2)x(N+2)
%%
%% The function works for a standard 1D domain, cylindrical
%% symmetric, or spherical symmetry depending on value set for
%% dim in the calling program.  
%% 
%% dim=0 for regular 1d
%% dim=1 for cylindrical symmetry
%% dim=2 for spherical   symmetry  
%%  
%% define matrices ML, MR
%%
  end1 = N+2;
  end2 = end1 + end1 - 1;
%% Diagonal
    j=1
    A(j,1)=j;
    A(j,2)=j;
    A(j,3)=-2*(1+dim);
    for j=2:end1
        A(j,1)=j;
        A(j,2)=j;
        A(j,3)=-( (x(j)+dx/2)^dim + (x(j)-dx/2)^dim  )/x(j)^dim;
    end;

%% Sub-Diagonal
    for j=2:end1
        A(end1+j-1,1)=j;
        A(end1+j-1,2)=j-1;
        A(end1+j-1,3)= ( (x(j)-dx/2)^dim  )/x(j)^dim;
    end;
%% Super-Diagonal
    for j=1:end1-1
        A(end2+j,1)=j;
        A(end2+j,2)=j+1;
        A(end2+j,3)= ( (x(j)+dx/2)^dim  )/x(j)^dim;
    end;
%%
%% For derivative BC at x=0 modify Super-Diagonal element in first row;
%%
    A(end2+1,3)=2*(1+dim);
%%    
%% For derivative BC at x=L modify Sub-Diagonal element in last row;
%%
    A(end1+end1-1,3)=A(end1+end1-1,3) + ( (x(end1)+dx/2)^dim  )/x(end1)^dim;
%%
    SA = spconvert(A)
    ML = speye(end1) - .5*lambda*SA;
    MR = speye(end1) + .5*lambda*SA;
return;