Solving Nonlinear Boundary Value Problem with Constraints

Question

I want to solve the following boundary value problem

where S, k_i, x_f, α_1, and θ are known parameters. We are trying to solve for h(x), p, and θ_d.

My idea was to use finite differences to create a numerical method for the problem. If we use second-order central differencing on the ODE, then we have the following equation:

Now the vanishing Dirichlet BCs tell me that h_0 = h_N = 0. However, I also need to use the contact angle boundary conditions to generate equations for p and \theta_d. To do this, I used a first-order central difference on the remaining BCs which gave me the following ghost point equations:

With these additional equations, I now should be able to apply the original scheme at all points i=0,...,N with the only caveat that h_0=h_N=0. Note that in this problem it is necessary to still include the equation of the scheme at i=0,N since these give us the additional information required to solve for the unknown constants p and θ_d.

Now my actual question is about the numerical implementation which I want to do in MATLAB since it is what I am most comfortable working with. I want to implement an iterative scheme that uses the old guess to update itself. As an initial guess, I will use the solution for the case of S=0 which can be solved for analytically. So at each iteration, I will need to solve a system of nonlinear equations. My first thought of how to do this was to use fsolve however this requires defining a function handle for each h_i. Is there any easy way of doing this?

As I was trying to think through the logic, I was able to figure out how this would work for the very simple case of N=4. With such a small number of points, one can actually just write out the equations into a function handle as I did in the code at the bottom. However, to accurately solve this, one of course needs way more points which makes writing out each equation very tedious. Is there a better way to do this? Is my logic even sensible?

clear; clc;
% Set problem parameters
xf = 6.17415;
theta = 39*(pi/180);
S = 1; ki=.37; alpha1 = 2.69;
% Set numerical parameters
L = 2*xf;
n = 4;
dx = L/n;
x=-xf:dx:xf;
iters = 5; % number of iterations

uold = zeros(1,n+3);   % First n-1 elements are h_1, ..., h_{n-1}; Last 2 elements are p and \theta_d
hold = (1 + (sinh(x-xf)-sinh(x+xf))/sinh(2*xf))*.809827;
uold(1,1:n+1) = hold;
uold(1,n+2:n+3) = [.809827,0];
plot(x,hold)
hold on
F = @(y) [1/(dx^2+3)*y(2)-(dx^2+2)/(dx^2+3)*y(1)-S/(dx^2+3)*exp(-2*ki*((x(2)+xf)+alpha1*hold(1,2)))+y(4)/(dx^2+3);
    1/(dx^2+3)*(y(3)+y(1))-(dx^2+2)/(dx^2+3)*y(2)-S/(dx^2+3)*exp(-2*ki*((x(3)+xf)+alpha1*hold(1,3)))+y(4)/(dx^2+3);
    1/(dx^2+3)*y(2)-(dx^2+2)/(dx^2+3)*y(3)-S/(dx^2+3)*exp(-2*ki*((x(4)+xf)+alpha1*hold(1,4)))+y(4)/(dx^2+3);
    y(4)+2*y(1)-dx*tan(theta-y(5))-S;
    y(5)+theta-atan(2/dx*y(3)-S/dx*exp(-4*ki*xf)+y(4)/dx);];

for j=1:iters
    u = fsolve(F,uold);
    plot(x,u(1,1:n+1))
    disp(u(1,n+2:n+3))
    uold = u;
end
hold off