# Numerical solution of PDE with uniform initial condition

I have a PDE like this

$$\frac{\partial h}{\partial t} = \bigg(\frac{\dot{L}}{L}\bigg)x\frac{\partial h}{\partial x} - \alpha\bigg[h^3\frac{\partial^3 h}{\partial x^3}\bigg]$$ With boundary and initial conditions given as

$$h(0,t) = h(1,t) = 1; h_{xxx}(0,t) = 0; h_{xxx}(1,t) = \frac{1}{\alpha}\frac{\dot{L}}{L}; h(x,0) = 1$$

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


Your PDE clearly states that if $$h$$ is uniform, it does not evolve (all spatial derivatives are zero). The only way to trigger an evolution is via the boundary conditions $$h_{xxx}(0,t)$$ and $$h_{xxx}(1,t)$$. It however seems that you never inject them in your duscretised PDE. Here your hxxx term at the boundaries does not involve the values of these imposed boundary conditions, since you only compute the gradients based in the inner points. My take would be to introduce ghost points (fake additional points, one on each side). Find the finite difference formula that uses them to compute $$h_{XXX}$$ on the sides, and use that to find what values of $$h$$ should be (dynamically) attributed to them so as to produce the desired values of $$h_{XXX}$$ at the sides.