function R = rhs(M,a,mu,dx,ap,am,u);

a1 = ap/dx + mu/dx^2;
a2 = -ap/dx + am/dx - 2.0*mu/dx^2;
a3 = -am/dx + mu/dx^2;

R(2:M) = a1*u(1:M-1) + a2*u(2:M) + a3*u(3:M+1) + 2*a;