How can I include a finite integral in a system of non-linear 2nd order differential equation using bvp4c function and solve numerically? Tried a test with Caughey example for large-amplitude whirling of an elastic string

In fine, I looking for a numerical solution using bvp4c of a 3rd order differential for a clamped circular membrane where the rhs term depends on a finite integral of the soolution y.

The test example is described in the book "Solving-odes-with-matlab" from Shampine Gladwell Thompson example 1.19 et 3.9, methodology is provide but no correction is given.

Quote ------

Caughy (1970) describes the large-amplitude whirling of an elastic string by a BVP consisting of the ODE:

with boundary conditions

Here α is a physical constant with 0 <α<1. Because the whirling frequency ω is to be determined as part of solving the BVP, there must be another boundary condition. Caughey specifies the amplitude ε of the solution at the origin: µ(0) = ε.

An unusual aspect of this problem is that an important constant H is defined in terms of the solution µ(x) throughout the interval of integration:

Formulate this BVP in standard form. As in the Sturm–Liouville example, you can introduce a new variable y3(x), a first-order ODE, and a boundary condition to deal with the integral term in the definition of H. The trick to dealing with H is to let it be a new variable y4(x). It is a constant, so this new variable satisfies the first-order differential equation y4' =0. It is given the correct constant value by the boundary condition resulting from the definition of H:

Unquote ------

I tried to code several ways, but always get one od the variable y(3) or y(4) out of bound. I miss some methodlogy to get a correct output. Thanks for any help.

The following code include following analysis for differential equation:

and

y4(x)=H=cste and

and boundary conditions

and

function dydx = bvpfcn(x,y,OMEGA)
  ALPHA=0.5 ; EPSILON=0.5;
  dydx = zeros(4,1);
  dydx = [y(2);...
          -OMEGA.^2.*(((1-ALPHA.^2)./y(4)).*(1./(1+y(1).^2).^0.5)+ALPHA.^2).*y(1);...
          (1./(1+y(1).^2).^0.5);...
          0];
end
function res = bcfcn(ya,yb,OMEGA)
  ALPHA=0.5 ; EPSILON=0.5;
  res = zeros(4,1);
  res = [ya(2);...
         yb(2);...
         yb(3)-(asinh(yb(1))-asinh(ya(1)));...
         yb(4)-((1./ALPHA.^2).*(1-(1-ALPHA.^2).*yb(3)))];
end
function g = guess(x)
  g = [cos(pi*x);-pi*sin(pi*x)];
end
xmesh = linspace(0,1,10);
OMEGA=pi; % initial condition
solinit = bvpinit(xmesh,@guess,OMEGA);
tol=1E-3;
options = bvpset("RelTol",tol,"AbsTol",tol,"Stats","on",'Nmax',10000); 
%
tic;
Testsol4c = bvp4c(@bvpfcn, @bcfcn, solinit);
toc;

[darova]

What is the question? You have a code. DOens't it work?

[Menfin]

Apreciate your reply Darova. The code does not work and return y(3) out of bound.