Adding Vector in discretization while solving Boundary Value Problem using bvp4c

I'm trying to solve a fourth order differential equation of the form:-

Here, p and W are both function of x, I calculate them from another function.

The 4 boundary conditions to solve this problem:- At x = 0, (1) (2)

At x = x(end) = L, (3) (4)

In this case, is a vector which is a function of x. At x=0, I say, , and at x = x(end),

Here's the code where I define the problem:-

options=bvpset('Abstol',2.56e-7,'Vectorized','on');

    function Y = guess(X)
        [~,E,I] = Sec_Prop(X,Section);
        Y = [0.003;-0.0015;mt/(E*I);pt/(E*I)-Px*0.0015/(E*I)];
    end

Solinit = bvpinit(x,@guess); %Initial condition for bvp4c.

S1= bvp4c(@odefun1,@bcfun1,Solinit,options); %Solution in the form of struct.

My code to define the ODE and Boundary Conditions:-

% Original Code
function dydx1 = odefun1(X,y1)
[D,E,I] = Sec_Prop(X,Section);
R = E.*I;
W = dist_load(X,W_Dist);
p = soilreact(X,D,y1(1,:));
dydx1 = [y1(2,:);y1(3,:);y1(4,:);(-p-Px.*y1(3,:)+W)./R];
end

function res1 = bcfun1(ya,yb)
res1 = [ya(4)+Px/(E_upper*I_upper)*ya(2)-pt/(E_upper*I_upper);ya(3)-mt/(E_upper*I_upper);yb(4)+Px/(E_lower*I_lower)*yb(2);yb(3)];
end

But I'm not getting the desired solution because the original Finite Difference equation to solve this problem looks like this:-

Where both y and R=E*I are being discretized. To my knowledge, only the y (the solution) is discretized by bvp4c.

I've somewhat solved the problem by rewriting my odefun and bcfun after assuming y1(1) = R*y, y1(2) = R*y' and so on. My modified code looks like this:-

% Modified Code
function dydx1 = odefun1(X,y1)
[D,E,I] = Sec_Prop(X,Section);
R = E.*I;
W = dist_load(X,W_Dist);
p = soilreact(X,D,y1(1,:)./R); % Here, y1(1,:)./R = y1(1,:) from the previos code, which is creating instability in my solution when R changes. 
dydx1 = [y1(2,:);y1(3,:);y1(4,:);-p-Px.*y1(3,:)./R+W];
end

function res1 = bcfun1(ya,yb)
res1 = [ya(4)+Px*ya(2)/(E_upper*I_upper)-pt;ya(3)-mt;yb(4)+Px*yb(2)/(E_lower*I_lower);yb(3)];
end

Which didn't solve the problem entirely. In this case, I'm having isssues on calculating p.

Can I get some idea to solve this issue, so that I can discretize both R and y in my solution, like the one done in the original solution? Thank You!