% Program 8.7 Implicit Newton solver for Burgers equation 
% input: space interval [xl,xr], time interval [tb,te], 
%        number of space steps M, number of time steps N 
% output: solution w 
% Example usage: w=p4_1(0,1,0,2,20,40) 
function w=p4_1(xl,xr,tb,te,M,N)
D=1;
f=@(x) .5+cos(2*pi*x);  
h=(xr-xl)/M; k=(te-tb)/N; m=M+1; n=N; 
sigma=D*k/(h*h); 
w(:,1)=f(xl+(0:M)*h)'; % initial conditions
w1=w; 
for j=1:n 
  for it=1:3 % Newton iteration 
    DF1=diag(2*k+1+2*sigma*ones(m,1))+diag(-sigma*ones(m-1,1),1);
    DF1=DF1+diag(-sigma*ones(m-1,1),-1);
    DF2=diag(-6*k*w1);
    DF3=diag(3*k*w1.^2);
    DF=DF1+DF2+DF3;
    F=-w(:,j)+(DF1+DF2/2+DF3/3)*w1;
    DF(1,:)=[-3 4 -1 zeros(1,m-3)];F(1)=DF(1,:)*w1;
    DF(m,:)=[zeros(1,m-3) -1 4 -3];F(m)=DF(m,:)*w1; 
    w1=w1-DF\F; 
  end 
  w(:,j+1)=w1; 
end 
x=xl+(0:M)*h;t=tb+(0:n)*k; 
mesh(x,t,w')
end