function [ C ] = Theta_C( theta,n,r ) %CN_C computes the matrix C such that u(j+1) = Cu(j) in the theta method % Input variables: % > theta: The value of theta for the theta method % > n: The size of the matrix % > r: The Crank-Nicolson constant % Output variable: % > C: The matrix C % First we require the base matrix that incorporates the boundary % conditions A = -2*eye(n); A = A + diag(ones(n-1,1),1) + diag(ones(n-1,1),-1); % Incorporate the Neumann boundary conditions A(1,1) = -1; A(n,n) = -1; % Use this to define the premultiplier matrix for u(j+1) AA = eye(n) - r*theta*A; % And define the premultiplier matrix for u(j) BB = eye(n) + r*(1-theta)*A; % Finally calculate C C = AA\BB; end