Numerical solution of PDE with uniform initial condition

I have a PDE like this

With boundary and initial conditions given as

In this case, I have a uniform initial condition and when I use the ode15s solver to solve this PDE, the final result I get is always 1. So, the plot looks like the picture below but the original results are supposed to be curvy as time increases. Not sure what I am doing wrong and I hope someone can help. I have attached the code below.

clear all
clc
%%
pts = 2^10-1; tmax = 1.76;
ti = linspace(0,tmax,pts);
xi = linspace(0,1,pts);
h_init = ones(pts,1);

[t,y] = ode15s(@(t,y) example_model(t,y,xi'),ti,h_init);

tspan = [1 100 200 300 400 500];
figure,
for i = tspan
    plot(xi,y(i,:));
    hold on;
end

function dhdt = example_model(t,y,xi)
t
pts = length(y); h = y;
h([1 pts]) = 1; h0 = 100;
dx = xi(2)-xi(1);
% L and Lt
tau = 0.765; U0 = 0.37;
lambda = 11.6; L_cl = 0.4;
tt = t./tau; st = sqrt(tt);
term1 = -0.5.*(tt).^2;
term2 = 0.5.*sqrt(pi).*erf(st);
term3 = -st.*exp(-tt);
L = L_cl + U0.*tau.*(term1 + lambda.*(term2+term3));
Lt = -U0.*tau.*(t/tau^2 - lambda.*((exp(-tt).*st)./tau - exp(-tt)./(2.*tau.*st) + (3991211251234741.*exp(-tt))./(4503599627370496.*tau.*pi^(1/2).*st)));
if isnan(Lt)
    Lt = 0;
end
% constant
C = 7.87e-7; alpha = C*h0^3/(3*L^4);
% derivatives
hx = FD_derivative(h,dx);
hxx = FD_derivative(hx,dx);
hxxx = FD_derivative(hxx,dx);
q = (h.^3).*hxxx;
qx = FD_derivative(q,dx);
dhdx = (Lt/L).*xi.*hx;
dhdt = dhdx - alpha.*qx;
end

function dydx = FD_derivative(y,dx)
N = length(y); dydx = zeros(N,1);
for i = 1:N
    if i == 1
        dydx(i) = -y(i) + y(i+1);
    elseif i == N
        dydx(i) = y(i) - y(i-1);
    else
        dydx(i) = -0.5*y(i-1) + 0.5*y(i+1);
    end
end
dydx = dydx./dx;
end

Answer 1

[VBBV]

clear all
clc
%%
pts = 2^0.1; tmax = 1.76;
ti = linspace(0,tmax,50);  % give suitable sample size 
xi = linspace(0,1,50);
h_init = ones(50,1);
[t,y] = ode15s(@(t,y) example_model(t,y,xi'),ti,h_init);
tspan = [1 100 200 300 400 500];
figure,
color = {'b','r','g','y','m','k'}

color = 1×6 cell array
    {'b'}    {'r'}    {'g'}    {'y'}    {'m'}    {'k'}

for i = 1:length(tspan)
    plot(xi,y(i,:),color{i});
    hold on;
end
legend('t = 1','t = 100','t =200','t = 300','t = 400','t = 500')

function dhdt = example_model(t,y,xi)
t;
pts = length(y); 
h = y;
h([1 pts]) = 0;  % change this boundary condition here 
h0 = 100;
dx = xi(2)-xi(1);
% L and Lt
tau = 0.765; U0 = 0.37;
lambda = 11.6; L_cl = 0.4;
tt = t./tau; st = sqrt(tt);
term1 = -0.5.*(tt).^2;
term2 = 0.5.*sqrt(pi).*erf(st);
term3 = -st.*exp(-tt);
L = L_cl + U0.*tau.*(term1 + lambda.*(term2+term3));
Lt = -U0.*tau.*(t/tau^2 - lambda.*((exp(-tt).*st)./tau - exp(-tt)./(2.*tau.*st) + (3991211251234741.*exp(-tt))./(4503599627370496.*tau.*pi^(1/2).*st)));
if isnan(Lt)
    Lt = 0;
end
% constant
C = 7.87e-7; alpha = C*h0^3/(3*L^4);
% derivatives
hx = FD_derivative(h,dx);
hxx = FD_derivative(hx,dx);
hxxx = FD_derivative(hxx,dx);
q = (h.^3).*hxxx;
qx = FD_derivative(q,dx);
dhdx = (Lt/L).*xi.*hx;
dhdt = dhdx - alpha.*qx;
end
function dydx = FD_derivative(y,dx)
N = length(y); dydx = zeros(N,1);
for i = 1:N
    if i == 1
        dydx(i) = -y(i) + y(i+1);
    elseif i == N
        dydx(i) = y(i) - y(i-1);
    else
        dydx(i) = -0.5*y(i-1) + 0.5*y(i+1);
    end
end
dydx = dydx./dx;
end

you can change the boundary condition to suit your model output inside the example model function. and give a reasonable step range at beginning

[Darsh]

Hi, but changing the boundary condition like you said changes the whole mathematical model for me. But I found my mistake. I forgot to introduce the third derivative boundary condition and it works better now. Thanks for the insight.

[VBBV]

Ok. I just showed that if you change the model inputs it will give more precise results.