% 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=ex13(0,1,0,1,10,10)
function w=ex13(xl,xr,tb,te,M,N)
%f=@(x) .5 + .5*cos(pi*x) % 8.13.a
f=@(x) .5 + .5*cos(pi*x)  %8.13.b
D = 1;
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(1-k+2*sigma*ones(m,1))+diag(-sigma*ones(m-1,1),1);
     DF1=DF1+diag(-sigma*ones(m-1,1),-1);
     DF2=diag(2*k*w1);
     DF=DF1+DF2;
     F=-w(:,j)+(DF1+DF2/2)*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