Matlab: Resolviendo la ecuación de onda bidimensional en el caso de una membrana de forma piramidal

El código fuente para resolver una ecuación de onda bidimensional está dado en el archivo twodimwavedirbc.m disponible en esta web: MATH 5343 Numerical Solution for Partial Differential Equations, Spring 2011.

La idea es resolver el ejercicio 12a del capítulo 12.2 del libro "Introduction to Numerical Ordinary and Partial Differential Equations Using MATLAB" de Alexander Stanoyevitch (de este libro proviene el scrip twodimwavedirbc.m):

Es decir, es una ecuación de onda cuya forma inicial es una pirámide. La dificultad aquí es dibujar la pirámide, yo la deduje como explico a continuacion:

La pirámide mostrada:

está formada por cuatro planos limitados por las rectas Y=X y Y = -X/2+4:

Las regiones donde Z es una función de X o de Y (dependiendo del plano donde se proyecta el eje Z) se muestran a continuación:

Para X<2, Y<2 y la región donde X<Y, la variable Z está definida por la función Z= Y/2.
Esto se hace para cada región de cada uno de los cuatro planos que definen la pirámide.

El script resultante es:

function z = pyramid(x, y)
z=0;

if y<=2 && x <= 2
   if y < x
       z=y/2;
   else
       z=x/2;
   end
end

if y<=2 && x > 2 && x <= 4
   if y < 4-x
       z=y/2;
   else
       z=-x/2+2;
   end
end

if y > 2 && x <= 2 && y <= 4
   if y < 4-x
       z=x/2;
   else
       z=-y/2+2;
   end
end

if y > 2 && y <= 4 && x > 2 && x <= 4
   if y < x
       z=-x/2+2;
   else
       z=-y/2+2;
   end
end

end 

Modifiqué un poco el archivo twodimwavedirbc.m para hacerlo más flexible:

function [x, y, t, U] = twodimwavedirbc2(phi, nu, a1, a, b1, b, T, h, c)
% a1: punto inicial del eje x
% b1: punto inicial del eje y

% solves the two-dimensional wave problem u_tt = c^2(u_xx+u_yy)
% on the rectangle {a1<= x <= a,  b1<= y <= b}, with u(x,y)=0
% on the boundary of the rectangle.

% variables de entrada: 
% phi=phi(x,y) = función de la forma inicial de la onda
% nu=nu(x,y) = función inicial de la velocidad onda

% a= right endpoint of x, b = right endpoint of y, 
% T = final time solution will be computed.  h=common gap on
% x, y-grids,   k= t-grid gap, c = speed of wave.
% wave.
% Output variables: x = row vector for first space variable, y =
% row vector for second space variable,  t = time grid row vector
%(starts at t=0, ends at t=T, has Nt equally spaced values), 
% U = (Nx)by(Ny)by(Nt) matrix of solution approximations at
% corresponding grid points (where Ny = number of y-grid points)
% x grid will correspond to first entries of U, y
% grid values to second entries of U, and t grid to third entries of U.
% CAUTION:  This method will only work if h is chosen so that the x
% and y grids can have a common gap size, i.e., if h = a/h and
% b/h must be integers.
% The time grid gap is chosen so that mu = 1/2, this guarantees the
% Courant-Friedrichs-Levy condition holds and simplifies the
% main finite difference formula.

k=h/sqrt(2)/c; %k is determined from mu = 1/2

if (b/h-floor(round(b/h))>10*eps)||(a/h-floor(round(a/h))>10*eps)
    fprintf('Space grid gap h must divide evenly into both a and b \r')
    fprintf(' Either try another input or modify the algorithm')
    error('M-file will exit')
end

Nx = round((a-a1)/h)+1; %number of points on x-grid
Ny = round((b-b1)/h)+1; %number of points on y-grid
Nt = floor(T/k)+1; %number of points on t-grid
U=zeros(Nx, Ny, Nt);
x=a1:h:a;
y=b1:h:b;
t=0:k:T;
% Recall matrix indices must start at 1.  Thus the indices of the
% matrix will always be one more than the corresponding indices that
% were used in theoretical development.

% Note that by default, zero boundary values have been assigned to all
% grid points on the edges of the rectangle (and, for the time being,
% at all grid points).
%Assign initial time t=0 values and next step t=k values.
for i=2:(Nx-1)
     for j=2:(Ny-1)
    U(i,j,1)=feval(phi,x(i),y(j));
    U(i,j,2)=.25*(feval(phi,x(i-1),y(j))+...
      feval(phi,x(i+1),y(j))+feval(phi,x(i),y(j-1))+...
      feval(phi,x(i),y(j+1)))+ k*feval(nu,x(i),y(j));     
     end
end
%Assign values at interior grid points
for ell=3:Nt %letter ell looks too much like number one
    U(2:(Nx-1),2:(Ny-1), ell) = ...
     +.5*(U(3:Nx,2:(Ny-1), ell-1)+ U(1:(Nx-2),2:(Ny-1), ell-1)...
     + U(2:(Nx-1),3:Ny, ell-1) + U(2:(Nx-1),1:(Ny-2), ell-1))...
     - U(2:(Nx-1),2:(Ny-1), ell-2);
end

Y el script principal es:

clc; clear all; close all;
phi = @(x, y)sin(2*x) *sin(2*y);
nu=@(x, y)0;
[x, y, t, U] = twodimwavedirbc2('pyramid', nu, 0, 4, 0, 4, 4, 0.1, 1);
[a b c] = size(U);
%ans = 26 26 46
for ell=1:c
     surf(x,y,U(:, :,ell));
     axis([0, 4, 0, 4 -1, 1]);
     M(:,:,ell) = getframe;
end
movie(M, 2, 4)


Lo he probado en Matlab, y acá está corriendo lo mejor que puede en Octave 4.4.1 a pesar de que me dice que la función Movie no está implementada: