% m43_BackwardEuler.m - order of accuracy study for backward
% Euler solution of heat eq.  u_t = u_xx, u(-1)=u(1)=0 

% Compute solution for various grid sizes:
  figure(1)
  kk = []; uu = [];
  for logh = 1:8
    h = 2^(-logh);
    x = (-1+h:h:1-h)';                 
    k = .25*h;  
    u = exp(-10*x.^4./(1-x.^2));       % smooth initial data
    L = length(x);
    % sparse matrix for implicit terms:
      a = (1+2*k/h^2); b = -k/h^2;
      main = a*sparse(ones(L,1));
      off  = b*sparse(ones(L-1,1));
      B = diag(main) + diag(off,1) + diag(off,-1);
    tmax = 1/8;
    nsteps = tmax/k;
    hold off, shg
    plt = plot(x,u,'linewidth',2);
    title('t = 0','fontsize',18)
    axis([-1 1 -.01 1.01]), grid, pause
    for step = 1:nsteps
      u = B\u;                         % backward Euler
      set(plt,'ydata',u)
      drawnow
    end
    kk = [kk; k]; uu = [uu; u(x==1/2)];
    title('t = 1/8','fontsize',18)
    pause
  end
  
% Plot deviations of u(x=1/2,t=1/8) from finest-grid value:
  figure(2), hold off
  I = 1:length(kk)-2;
  loglog(kk(I),kk(I),'--r')
  text(.013,.025,'k','fontsize',18,'color','r')
  axis([.003 .15 3e-6 .5]), grid on, hold on
  loglog(kk(I),kk(I).^2,'--r'), grid on
  text(.013,.0004,'k^2','fontsize',18,'color','r')
  loglog(kk(I),abs(uu(I)-uu(end)),'.-k','markersize',18)
  xlabel('time step size k','fontsize',16)
  ylabel('error at  t=0.1','fontsize',16)