% m40_fourthorderdiffusion.m - explicit 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, 4e-8, 5e-8, 4.88e-8, 4.90e-8) '); 

% Sparse matrix for finite differencing:
  L = length(x);
  a = (1-6*k/h^4); 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));
  A = diag(main) + diag(off,1) + diag(off,-1) ...
                 + diag(off2,2) + diag(off2,-2);
  
% Time-stepping:
  t = 0;
  while 1
    u = A*u; t = t+k; set(plt,'ydata',u)
    title(['t = ' num2str(t)],'fontsize',20)
    drawnow
  end