function [D0,D2,D4,Mt,Dl,Dr]=init(nx,dx,dt,mu,alpha,c1)

e=ones(nx,1);

D0=spdiags([e],0,nx,nx);
D2=1/(dx*dx)*(spdiags([e -2*e e],-1:1,nx,nx)+spdiags([e],nx-1,nx,nx)+spdiags([e],-nx+1,nx,nx));
D4=1/(dx*dx*dx*dx)*(spdiags([e -4*e 6*e -4*e e],-2:2,nx,nx)+spdiags([e -4*e],nx-2:nx-1,nx,nx)+spdiags([-4*e e],-nx+1:-nx+2,nx,nx));
Mt=1/dt*spdiags([e],0,nx,nx);

Dl=Mt-alpha*(1-i*c1)*D2-alpha*(mu-1)*D0;
Dr=Mt+(1-alpha)*(1-i*c1)*D2+(1-alpha)*(mu-1)*D0;

