Resolviendo una ecuación parabólica unidimensional en Matlab/Octave

El siguiente ejercicio se encuentra en el libro "Computational Partial Differential Equations Using MATLAB" by Yi-Tung Chen, Jichun Li, capítulo 2, ejercicio 3:

Una pista de cómo resolver este ejercicio está en la sección 2.21:

Los autores dan un código base sobre el qué trabajar. Luego de recibir un poco de ayuda de la gentita de Stackoverflow, este es el código para resolver el problema:


 %---------------------------------------------------------------
% para1d.m:
%    use the explicit scheme to solve the parabolic equation
%    u_t(x,t) = u_{xx}(x,t),          xl < x < xr, 0 < t < tf
%    u(x,0) = f(x),                   xl < x < xr
%    u_x(0,t) = gl(t), u_x(1,t) = gr(t),  0  < t < tf
%---------------------------------------------------------------
clc;
close all;
clear all;                  % clear all variables in memory

xl=0; xr=1;                 % x domain [xl,xr]
J = 31;                     % J: number of division for x
dx = (xr-xl) / J;           % dx: mesh size
tf = 0.1;                   % final simulation time
Nt = 51;                    % Nt: number of time steps
dt = tf/Nt/4;

mu = dt/(dx)^2;

if mu > 0.5         % make sure dt satisy stability condition
    error('mu should < 0.5!')
end

% Evaluate the initial conditions
x = xl : dx : xr;              % generate the grid point
f = x;

% store the solution at all grid points for all time steps
u = zeros(J+1,Nt);
u_ex = zeros(J+1,Nt);

% boundary condition at left side
gl = 0;
% boundary condition at right side
gr = 0;

% Find the approximate solution at each time step
for n = 1:Nt
    t = n*dt;         % current time
     
    if n==1
      for j=2:J
       % first time step
        u(j,1) = f(j) + mu*(f(j+1)-2*f(j)+f(j-1));
      end
    end
  
    if n > 1
        for j=2:J
         % interior nodes
           u(j,n) = u(j,n-1) + mu*(u(j+1,n-1) - 2*u(j,n-1) + u(j-1,n-1));
        end
    end
  
    % hint 2.21 Los valores en la frontera es mejor calcularlos al final
    % pues estoy cogiendo índices que dependen de los nodos interiores
    u(1, n) = u(2, n) - dx * gl;
    u(J+1,n) = u(J, n) + dx * gr;
  
    % calculate the analytic solution
    for j=1:J+1
      xj = xl + (j-1)*dx;    
      suma = zeros(1000 , 1);

      for k= 1:1000
        suma(k) = 4/(((2*k-1)*pi)^2);
        suma(k) = suma(k) * exp(-((2*k-1)^2)*t*pi^2) * cos((2*k-1)*pi*xj);
      end

      u_ex(j, n)= double(0.5 - sum(suma));
    end
end

% Plot the results
tt = dt : dt : Nt*dt;
figure(1)
colormap(cool);     % draw gray figure
surf(x, tt, u');     % 3-D surface plot
xlabel('x')
ylabel('t')
zlabel('u')
title('Numerical solution of 1-D parabolic equation')

figure(2)
surf(x, tt, u_ex');     % 3-D surface plot
xlabel('x')
ylabel('t')
zlabel('u')
title('Analytic solution of 1-D parabolic equation')

maxerr=max(max(abs(u-u_ex))),


El resultado es:

Probado en Octave 4.4.1 y Matlab 2008 R2.