function w=Burgers(xl,xr,tb,te,M,N)
% Ecuación de Burges no lineal mediante por diferencias atrasadas
%   Entrada: [xl,xr], tiempo[tb,te], pasos M, tiempos N
 alfa=5; beta=4; D=0.05;
 f=@(x) 2*D*beta*pi*sin(pi*x)./(alfa+beta*cos(pi*x));
 l=@(t) 0*t; r=@(t) 0*t; 
 m=M+1; n=N; close all
 h=(xr-xl)/M; k=(te-tb)/N; sigma=D*k/h/h;
 w(:,1)=f(xl+(0:M)*h)'; w1=w;
 for j=1:n                   
    for it=1:4
       DF1=diag(1+2*sigma*ones(m,1))+diag(-sigma*ones(m-1,1),1);
       DF1=DF1+diag(-sigma*ones(m-1,1),-1);
       DF2=diag([0;k*w1(2:m-1)/(2*h)],1)-diag([k*w1(2:m-1)/(2*h);0],-1);
       DF=DF1+DF2;
       F=-w(:,j)+(DF1+DF2/2)*w1;
       DF(1,:)=[1 zeros(1,m-1)]; DF(m,:)=[zeros(1,m-1) 1];
       F(1)=w1(1)-l(j); F(m)=w1(m)-r(j);
       w1=w1-DF\F;                      % Newton-Raphson
    end
    w(:,j+1)=w1;
 end
 x=xl+(0:M)*h; t=tb+(0:n)*k; mesh(x,t,w')
end