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

R = rhs(M,a,mu,dx,ap,am,u);	
  
resid = norm(R);
  
a1 = -ap/dx - mu/dx^2;
a2 = +ap/dx - am/dx + 2.0*mu/dx^2;
a3 = am/dx - mu/dx^2;

a(2:M) = h*a1;
b(2:M) = 1 + h*a2;
c(2:M) = h*a3;

du = 0.0*u;
du = trib(a,b,c,h*R,2,M);
du(M+1) = 0.0;

u = update(M,1.0,u,du);