function [ML,MR,R] = matrixdefine(N,lambda)
%%
%% (N+1)*dx = L
%% For dirichlet at x=0 and neumann at x=L
%%     xj = (j-1)*dx so j=1 corresponds to x=0
%%                  and j=N+2 corresponds to x=L.
%% There are N+1 unknowns (at j=2:N+2), so matrices
%%     need to be (N+1)x(N+1)
%%  
%% define matrices ML, MR
%%
%% Diagonal
    for j=1:N+1
        A(j,1)=j;
        A(j,2)=j;
        A(j,3)=-2;
    end;

%% Sub-Diagonal
    for j=2:N+1
        A(j+N+1-1,1)=j;
        A(j+N+1-1,2)=j-1;
        A(j+N+1-1,3)=1;
    end;

%% Super-Diagonal
    for j=1:N+1-1
        A(j+2*(N+1)-1,1)=j;
        A(j+2*(N+1)-1,2)=j+1;
        A(j+2*(N+1)-1,3)=1;
    end;

%%C = -2*eye(N+1) + diag(ones(N,1),1) + diag(ones(N,1),-1);

%% First define for Neumann then Dirichlet conditions
%% For Neumann   at x=L set   A(2*(N+1)-1,3)=2;
    A(2*(N+1)-1,3)=2;
%% For Dirichlet at x=0 leave A(2*(N+1),3) unchanged;
    SA = spconvert(A);
    ML = speye(N+1) - .5*lambda*SA;
    MR = speye(N+1) + .5*lambda*SA;
%% For Neumann at x=L divide last row or ML and MR by 2 
%% so that ML will be symmetric for PCG
       ML(N+1,N  ) = ML(N+1,N  )/2.;
       ML(N+1,N+1) = ML(N+1,N+1)/2.;
       MR(N+1,N  ) = MR(N+1,N  )/2.;
       MR(N+1,N+1) = MR(N+1,N+1)/2.;
    R = cholinc(ML,1e-3);
return;