pdepe for a fourth-order de

Hi everyone. I'm trying to solve the following pde:

and boundary conditions:

Basically, this is a "gradient flow" method for solving classical ODEs; in this case I want to solve the equation u_xxxx = f(x). Henceforth, I build by hand the term F(x) whenever the solution u(x) is chosen.

global F
%%MESH
xa = 0;
xb = 1;
nx = 100;
xmesh = linspace(xa,xb,nx);

%%FUNCTIONS
syms x
U0 = 2*x^2 - x^4 + exp(x)/10 + sin(5*pi*x)/10;
U1 = diff(U0,1);
U2 = diff(U0,2);
U3 = diff(U0,3);
U4 = diff(U0,4);

U0 = matlabFunction(U0);
U1 = matlabFunction(U1);
U2 = matlabFunction(U2);
U3 = matlabFunction(U3);
U4 = matlabFunction(U4);

F = U4;

Using pdepde, I set u'' = y in order to get a lower-order pde, then:

Up to now, this is my code:

coefficients

function [c,f,s] = funpde1coeff(x,t,u,DuDx)
global F
%
c = [1; 0]; 
% 
f = [0 -1;1 0] * DuDx; 
%
s = [F(x);-u(2)];

initial conditions

function u0 = funpde1ic(x)
global G
u0 = [G(x); 0];

where the function G(x) is built in such a way it satisfies the BCs (basically it is a third order polynomial function with 4 parameters).

Now how can I impose the original boundary conditions on the first derivative? Moreover, what if I have BCs on both first and second derivative at edges?

Thanks in advance for your help. Best,

Pinco

%============= EDIT ===========

I tried with BCs on second derivative, but I get the following error:

Warning: Failure at t=0.000000e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (7.905050e-323) at time t.

The code is:

global U0a U1a U2a U0b U1b U2b 
U0a = U0(xa);
U1a = U1(xa);
U2a = U2(xa);
U0b = U0(xb);
U1b = U1(xb);
U2b = U2(xb);

global F2
F2 = U4;
tspan2 = linspace(0,0.2,10);
m2 = 0;
SOL2 = pdepe(m2,@funpde2coeff,@funpde2ic,@funpde2bc,xmesh,tspan2);

function [c,f,s] = funpde2coeff(x,t,u,DuDx)
global F2
%
c = [1; 0]; 
% 
f = [0 -1;1 0] * DuDx; 
%
s = [F2(x);-u(2)]; 

function u0 = funpde2ic(x)
global G2
u0 = [G2(x); 0];

function [pl,ql,pr,qr] = funpde2bc(xl,ul,xr,ur,t)
global U0a U0b U2a U2b
% left
pl = [ul(1) - U0a ; ul(2)- U2a]; 
ql = [0; 0]; 
% right
pr = [ur(1) - U0b ; ur(2)- U2b];
qr = [0; 0];

Where is the mistake?

Answer 1

[Pinco]

any ideas?

[Torsten]

Why don't you use bvp4c to solve u_xxxx=f(x) directly ?

Best wishes

Torsten.