function [ctempnew] = matsolve(ML,MR,ctempold,bcvector);
%%
%% Solve ML*ctempnew = MR*ctempold + (lambda/2) c(x=0,old) + (lambda/2) c(x=0,new).
%% bcvector holds (lambda/2) c(x=0,old) + (lambda/2) c(x=0,new).
%% Use preconditioned conjugate gradient with incomplete Choleski
%% preconditioning to solve the linear system.  initial guess: c_old
RHS = MR*ctempold + bcvector;
ctempnew=ML\RHS;
return;