Using array-defined coefficient in bvp4c

Hello Matlab Community,

currently I try to solve a boundary value problem using bvp4c. If the coefficients in the ODE are constant or described analytical as function of x (here the first coefficient "-phi2*x*y(1)"), the solution works:

solinit = bvpinit(linspace(0,lBeam,10),[0 0 0 0]);
sol     = bvp4c(@twoode, @twobc, solinit);
y_num   = deval(sol, x_num);

with

function [dydx] = twoode2(x, y)

      global mu2 phi2 alpha_l alpha_l2

      dydx = [y(2); y(3); y(4); -phi2*x*y(1)+mu2*y(3)+mu2*alpha_l2*cos(alpha_l*x)];
  end

and

function [res] = twobc(ya,yb)
% Boundary Condition

      res=[ya(1); yb(1); ya(3); yb(3)];

end

But for a better approximation of reality I need to put in a parameter "D(x)" which is not representable as analytical function but as section-defined function (array). So dydx change to the following:

   dydx = [y(2); y(3); y(4); -phi2*D*y(1)+mu2*y(3)+mu2*alpha_l2*cos(alpha_l*x)];

with D for example (length is not interesting): start analysis

D=[1 1 1 1 1 1 1 1...]; % length depends on the points used to solve the BVP

adapting D after analysis

D=[1 0.1 1 1 0.1 1 0.1 1...]; % value 0.1 depends on position x

So there is no function to describe D because it depends on the solution of the ODE and can be adapted from solution to solution.

The question is not how to adapt D but to put in an array-coefficient in the ODE. If D=x (linear function) the algorithm bvp4c also just uses points (here array x) to get solution.

I hope someone can help me!

Flori