ode15s solver produces wrong solution

I am trying to solve the heat equation

such that

using the method of line. I discritzed in space and used ode15s to integrate the semi-discrete system of ODEs using ode15s in time, but still when I compare to the exact solution it is wrong. Is there anything wrong in my code?

t0 = 0; % initial time
Tf = 1; % final time
xl = 0; % left point
xr = 1; % right point
m = 52; % no. of ODEs in semidiscrete system
h = 1/m; % space stepsize
k = 0.0001; % time stepsize
Utrue = @(t,x) exp(-t/10) * sin(pi*x); % true solution
xsol = linspace(xl,xr,xr/h);
tspan = linspace(t0,Tf,Tf/k+2);
% initial conditions
U0=zeros(m,1);
for i=1:m
    x=i*h;
    U0(i) = Utrue(0,x);
end
% solve
opts = odeset('RelTol',1e-13, 'AbsTol',1e-13*ones(m,1),'InitialStep',k, 'MaxStep',k);
[t,u]=ode15s(@ODE,tspan,U0,opts);

% plot solution at x = 0.5 >>> Wrong Solution!!
figure(1)
plot(t,u(:,m/2))
hold on
utrue = Utrue(t,1/2);
plot(t,utrue)

% supporting function
function Ut = ODE(t,U)
m = 52;
h = 1/m; % space stepsize
% BCs
g = zeros(m,1); 
Utrue = @(t,x) exp(-t/10) * sin(pi*x); % true solution u
G = @(t,x)  -0.1*exp(-t/10)*sin(pi*x) + exp(-t/10)*pi^2 * sin(pi*x); % G = ut-uxx
g(1) = Utrue(t,0)/h^2 + G(t,0);
g(m) = Utrue(t,1)/h^2 + G(t,1);
%semidiscrete system
Ut = zeros(m,1);
K = 1;
Ut(1) = (K/h^2)*(-2*U(1) + 1*U(2));
for i = 2:m-1
    Ut(i) = (K/h^2)*(U(i-1) - 2*U(i) + U(i+1));
end
Ut(m) = (K/h^2)*(U(m-1) - 2*U(m));
Ut = Ut + g;
end

Answer 1

[Torsten]

t0 = 0; % initial time
tf = 1; % final time
nt = 20; % number of output times
xl = 0; % left point
xr = 1; % right point
nx = 51; % no. of ODEs in semidiscrete system
dx = 1/(nx-1); % space stepsize
Utrue = @(t,x) exp(-t/10).*sin(pi*x); % true solution
x = linspace(xl,xr,nx);
tspan = linspace(t0,tf,nt);
% initial conditions
U0 = Utrue(0,x);
% solve
%opts = odeset('RelTol',1e-13, 'AbsTol',1e-14*ones(m,1));
[t,u]=ode15s(@(t,y)ODE(t,y,nx,dx,x),tspan,U0);%,opts);
% plot solution at x = 0.5 
figure(1)
plot(t,u(:,(nx+1)/2))
hold on
utrue = Utrue(t,1/2);
plot(t,utrue)

% supporting function
function Ut = ODE(t,U,nx,dx,x)
  G = @(t,x)  -0.1*exp(-t/10).*sin(pi*x) + exp(-t/10).*pi^2.*sin(pi*x); % G = ut-uxx
  %semidiscrete system
  K = 1;
  Ut = zeros(nx,1);
  %Ut(1) = K*(-2*U(1) + 1*U(2))/h^2 + G(t,0);
  Ut(2:nx-1) = K*(U(1:nx-2) - 2*U(2:nx-1) + U(3:nx)) / dx^2 + G(t,(x(2:nx-1)).');
  %Ut(nx) = K*(U(nx-1) - 2*U(nx))/dx^2 + G(t,1);
end

[Romio]

@Torsten Thank you very much. I am not sure though how the boundary conditions at x = 0 and x =1 are treated? I know for this problem it is u(t,0) = u(t,1) = 0, but what if it is not?

[Torsten]

I know for this problem it is u(t,0) = u(t,1) = 0, but what if it is not?

Then you have a different Utrue (and a different G).

If the boundary conditions at x=0 and x=1 would change with time or whatever, your function "ODE" would be

function Ut = ODE(t,U,nx,dx,x)
  Utrue = @(t,x) exp(-t/10).*sin(pi*x); % true solution
  G = @(t,x)  -0.1*exp(-t/10).*sin(pi*x) + exp(-t/10).*pi^2.*sin(pi*x); % G = ut-uxx
  %semidiscrete system
  U(1) = Utrue(x(1),t);
  U(nx) = Utrue(x(nx),t);
  K = 1;
  Ut = zeros(nx,1);
  %Ut(1) = K*(-2*U(1) + 1*U(2))/h^2 + G(t,0);
  Ut(2:nx-1) = K*(U(1:nx-2) - 2*U(2:nx-1) + U(3:nx)) / dx^2 + G(t,(x(2:nx-1)).');
  %Ut(nx) = K*(U(nx-1) - 2*U(nx))/dx^2 + G(t,1);
end

Of course, this function also works in your present case, but is not necessary because Ut(1) = Ut(nx) = 0 is set and U(1) = U(nx) = 0 is already coming from the initial conditions for U.