Implementing numerical method for PDE

Hello

I am trying to solve the following PDE with intital and boundary conditions such that . I used the second order centered finite difference discrtization for and then want solve the ode system using ode15s

Here is my attempt. When I plot the solution obtained from ode15s and compare it to the exact solution they are different. I am not if I made a mistake somewhere. Help is really appreciated

clc,clear,close all

% parameters
t0 = 0;
T = 1.0;
tspan = [t0 T];
xl = 0;
xr = 1;
m = 20;
x = linspace(xl,xr,m + 1);
dx = 1/m;
Uexact = @(t,x) exp(1i*(x-t));

% initial conditions
U0 = Uexact(0,x)';
U0 = U0(2:end-1);

% solve
fn = @(t,U) ODE(t,U,m,dx);
opts = odeset('RelTol',1e-13, 'AbsTol',1e-15);
[t,U] = ode15s(fn, tspan, U0, opts);

%compare with exact solution
plot(x(2:end-1),U(end,:))
hold on 
plot(x(2:end-1),Uexact(T,x(2:end-1)))


 
function dUdt = ODE(t,U,m,dx)
    A = eye(m-1);
    A = A * (-2);
    A = A + diag(ones(m-2,1),1);
    A = A + diag(ones(m-2,1),-1);
    A = (1/dx^2) * A;
    g = zeros(m-1,1);
    g(1) =   g(1)   + (1/dx^2) * exp(1i*(-1*t));
    g(end) = g(end) + (1/dx^2) * exp(1i*(1-t));
    dUdt = (1i) * (A*U) + g;
end

Thanks

[Davide Masiello]

What are exactly the boundary conditions?

[Davide Masiello]

Nevermind, I got it, see answer below.

Answer 1

[Davide Masiello]

clear,clc
tspan   = [0 1];
N       = 100; 
x       = linspace(0,1,N);
dx      = 1/(N-1);
Uexact  = @(t,x) exp(1i*(x-t));

U0      = Uexact(0,x);

M       = eye(N);
M(1,1)  = 0;
M(N,N)  = 0;
opts    = odeset('Mass',M,'MassSingular','yes');

[t,U]   = ode15s(@(t,U)yourPDE(t,U,N,dx), tspan, U0, opts);

plot(x,real(U(end,:)),'k',x(1:4:end),real(Uexact(1,x(1:4:end))),'r.')
xlabel('x')
ylabel('U')
title('At final time')
legend('Numerical','Exact','Location','Best')

plot(x,real(U(1:3:end,:)),'k',x(1:3:end),real(Uexact(t(1:3:end),x(1:3:end))),'r.')
xlabel('x')
ylabel('U')
title('At several times')

  
function dUdt = yourPDE(t,U,N,dx)
dUdt(1,1)     = U(1)-exp(-1i*t);    
dUdt(2:N-1,1) = 1i*(U(1:end-2)-2*U(2:end-1)+U(3:end))/dx^2;
dUdt(N,1)     = U(end)-exp(1i*(1-t)); 
end