function [u,resid] = PC(M,a,mu,dx,ap,am,h,u)

R = rhs(M,a,mu,dx,ap,am,u);	
  
resid = norm(R);
  
up = update(M,0.5*h,u,R);

R = rhs(M,a,mu,dx,ap,am,up);	
  
u = update(M,h,u,R);