Why is ode45 failing to solve this system of ODEs?

Using the Chebyshev's differentiation matrix to do spatial derivatives and the ode45 function for solving temporal ODE's, I am trying to numerically solve the following heat equation:

with the initial condition and boundary conditions . Note that and are Chebishev-Lobatto points and is the first root of the Bessel's function of first kind, order zero. To do this, I use the following cheb function from the book Spectral Methods in Matlab written by Lloyd N. Trefethen (not formatted to keep it as an exact copy from the book):

% CHEB compute D = differentiation matrix, x = Chebyshev grid
function [D,x] = cheb(N)
if N==0, D=0; x=1; return, end
x = cos(pi*(0:N)/N)';
c = [2; ones(N-1,1); 2].*(-1).^(0:N)';
X = repmat(x,1,N+1);
dX = X-X';
D = (c*(1./c)')./(dX+(eye(N+1))); % off-diagonal entries
D = D - diag(sum(D')); % diagonal entries
end

and I wrote the following code:

close all
clear
clc

r_max = 100;
N     = 20;

[D, x] = cheb(N);

y = x';
[yy, xx] = meshgrid(y, x);

% 2.4048 is the first root of the Bessel function of first kind, order zero
u_0 = besselj(0, 2.4048 * x) * sin(pi * y);
u_0 = u_0(:);

xx = xx(:);
yy = yy(:);

D([1, end], :) = 0; 
D(:, [1, end]) = 0; 

I  = eye(N + 1);

Dx  = kron(I, D);
Dxx = kron(I, D^2);
Dyy = kron(D^2, I);

L = (Dxx + Dx ./ (xx + 1)) ./ (r_max^2) + Dyy;

[t, u] = ode45(@(t, u) odefcn(t, u, L), [0, 1.5], u_0);

uu = zeros(length(t), N + 1, N + 1);

for k = 1 : length(t)
    
    uu(k, :, :) = reshape(u(k, :), N + 1, N + 1);
end

t_step = 2;

figure('Name', 'Solution to the heat equation'); 
surface(y, x, squeeze( uu(t_step, :, :) ), 'edgecolor', 'none');
xlabel('y'); 
ylabel('x'); 
grid on
colormap(jet(256)); 
title('Heat function');
colorbar; 

where the odefcn I wrote is:

function u_t = odefcn(t, u, L)

    u_t = L * u;

end

Since I got a lot of NaN elements in my resulted heat function, I started debugging. I made sure my Laplacian was OK using this code:

close all
clear
clc

r_max = 30;
N     = 100;

[D, x] = cheb(N);

y = x';

u    = sin( 8 * (x - 1) ) * (y.^3);
u_x  = 8  * cos(8 - 8 * x) * (y.^3);
u_xx = 64 * sin(8 - 8 * x) * (y.^3);
u_yy = 6  * sin(8 * (x - 1) ) * y;

L_u_a = ( u_xx + u_x ./ (x + 1) ) ./ (r_max^2) + u_yy; % Analitic Laplacian

[yy, xx] = meshgrid(y, x);

xx = xx(:);
yy = yy(:);

I  = eye(N + 1);

Dx  = kron(I, D);
Dxx = kron(I, D^2);
Dyy = kron(D^2, I);

L = (Dxx + Dx ./ (xx + 1)) ./ (r_max^2) + Dyy;

Lu = L * u(:);

L_u_n = reshape(Lu, N + 1, N + 1); % Numeric Laplacian

figure('Name', 'Analitic Laplacian'); 
surf(y, x, L_u_a, 'edgecolor', 'none');
xlabel('y'); 
ylabel('x'); 
grid on
colormap(jet(256)); 
title('Analitic Laplacian');
colorbar; 

figure('Name', 'Numeric Laplacian'); 
surf(y, x, L_u_n, 'edgecolor', 'none');
xlabel('y'); 
ylabel('x'); 
grid on
colormap(jet(256)); 
title('Numeric Laplacian');
colorbar; 

figure('Name', 'Error'); 
surface(y, x, L_u_a - L_u_n, 'edgecolor', 'none');
xlabel('y'); 
ylabel('x'); 
grid on
colormap(jet(256)); 
title('Error');
colorbar; 

which means I am making a mistake in the way I am using the ode45 function. I checked here if I am passing the arguments like I should. Seems to me like I am.

Can someone please explain to me why am I using the ode45 function incorrectly and how can I fix the arguments of it, as well as my odefcn file? This is actually the first time I am using it in this way (with an input argument L). Before I used it to solve systems of ODE's, but the numeber of those equations was usually up to four. Now, I am trying to solve a system with a lot of ODE's. Could I have exceedded the limit?

[John D'Errico]

Could you have exceeded the limit? No.

Answer 1

[Torsten]

Your matrix L has 9162 NaN values.

This is most probably caused by the division by (x+1) which is problematic for the problem being posed on [-1,1] x [-1,1].

Further, heat conduction problems are generally stiff. ode15s is the solver of choice in this case.

[Nikola Ristic]

Thank you for the answer. Seems to fairly work good after I removed the singularities and switched to ode15s. :D