% m41_implicit.m - implicit solution of 4th-order eq
%  u_t = -u_xxxx,   u(-1) = u(1) = u'(-1) = u'(1) = 0.

% Grid, initial data, and plotting setup:
  h = .025; x = (-1+h:h:1-h)';                 
  u = abs(x)<.3;             
  hold off, shg
  plt = plot(x,u,'linewidth',2);
  axis([-1 1 -.1 1.15]), grid
  k = input('k? (e.g. 1e-8, 1e-5, .01) ');

% Sparse matrix for finite differencing:
  L = length(x);
  a = (1+6*k/h^4);            %  note that the signs have changed
  b = -4*k/h^4; c = k/h^4;
  main = a*sparse(ones(L,1));
  off  = b*sparse(ones(L-1,1));
  off2 = c*sparse(ones(L-2,1));
  B = diag(main) + diag(off,1) + diag(off,-1) ...
                 + diag(off2,2) + diag(off2,-2);
  
% Time-stepping:
  t = 0;
  while 1
    u = B\u;              %  note B\u, not A*u
    t = t+k; set(plt,'ydata',u)
    title(['t = ' num2str(t)],'fontsize',20)
    drawnow
  end