% illustration of the heat equation
% Solve the heat equation using finite differences and Forward Euler
function main()
 
   % the number of data points. More points means prettier picture.
   N = 400;
 
   L = 2.5; % the box size is [-L, L] x [-L, L]
 
   XX = linspace(-L, L, N);
   YY = linspace(-L, L, N);
   [X, Y] = meshgrid(XX, YY);
 
   scale = 2;
   Z = get_step_function (N, scale, X, Y);
 
   CFL = 0.125; % CourantFriedrichsLewy
   dx = XX(2)-XX(1);  dy = dx; % space grid
   dt = CFL*dx^2;
 
   plot_dt = 0.004; % plot every plot_dt iterations
 
   cmap = rv_matrix_rows(autumn); % colormap
 
   % Solve the heat equation with zero boundary conditions
   T = 0:dt:1;
   iter = 0;
   frame_no = 0;
   for t=T
 
      % plot the current temperature distribution
      if floor(t/plot_dt) + 1 > frame_no
 
         frame_no = frame_no + 1
 
        % plot the surface
         figure(2); clf; 
         surf(X, Y, Z);
 
         %  make the surface beautiful
         shading interp; colormap(cmap); 
 
         % add in a source of light
         camlight (-50, 54);
         lighting phong;
 
         % viewing angle
         view(-40, 38);
 
         axis equal; axis off;
         axis([-L, L, -L, L, 0, scale])
 
         hold on; plot3(0, 0, 3.4, 'g*'); % a marker to help with cropping
 
         pause(0.1);
         %return
 
         file = sprintf('Movie_frame%d.png', 1000+frame_no);
         %saveas(gcf, file) %save the current frame
         print(gcf, '-dpng', '-r400', file) %save the current frame
 
         disp(file); %show the frame number we are at
 
         % cut at max_fr_no frames
         max_fr_no = 15; 
         if frame_no >= max_fr_no
            break
         end
 
      end
 
      % advance in time
      W = 0*Z;
      for i=2:(N-1)
         for j=2:(N-1)
 
            W(i, j) = Z(i, j) + dt * ( Z(i+1, j) + Z(i-1, j) + Z(i, j-1) + Z(i, j+1) - 4*Z(i, j))/dx^2;
 
         end
      end
      Z = W;
 
   end
 
 
% The gif image was creating with the command 
% convert -antialias -loop 10000  -delay 20 -compress LZW Movie_frame10* Heat_eqn.gif 
 
% get a function which is 1 on a set, and 0 outside of it
function Z = get_step_function(N, scale, X, Y)
 
   c = 2;
   d=-1;
   e=1;
   f=0.5;
   k=1.2;
   shift=10;
 
   Z = (c^2-(X/e-d).^2-(Y/f).^2).^2 + k*(c+d-X/e).^3-shift;
 
   Z = 1-max(sign(Z), 0);
   Z = scale*Z;
 
function X = rv_matrix_rows(X)
 
   [m, n] = size(X);
 
   for i = 1:m
      j = m + 1 - i;
      if i < j
         tmp = X(i, :); X(i, :) = X(j, :); X(j, :) = tmp;
      end
   end