How do I implement shooting method with ode solvers in Matlab?

Hi, I am trying to implement a shooting method with ode15s solver. I have a set of 4 non linear coupled differential equations. I know the initial value of three of the equations and I know the boundary condition for my problem. I have added images of the equations and initial conditions.

When I use fzero to guess the initial value of Q, the ode solver never goes till the end of the specified rspan.

I have also attached results from my simulation and the expected result. Any help with implementing a shooting method would be appreciated.

Thanks,

Sid

Code:

   function shooting_microfilm19

clear
close all

x0 = 0;

f = @(x) shooting_microfilm_thickness(x);
optionsf = optimset('display','off');
[xopt,fopt] = fzero(f,x0,optionsf);

disp(fopt)
disp(xopt)
end

function Krel=shooting_microfilm_thickness(x)%(M,R_g0,Tlv,rho_v,hlv,Tv,rho_l,sTen,Tsat,DT_sup,th0,nu_l,k_l,Tw,Pv,R_r,Rb)

% InputsWater1atm;
%water
%Properties at 100 deg C
Tsat = 373; %saturation temp in K
DT_sup = 5; %wall superheat
Tw = Tsat+DT_sup; %wall temp in K
% Vapour properties

rho_v=0.59814; % vapour density kg/m3

% Liquid properties
M=0.018; %water molecular wt
rho_l=958.35;%density of liquid kg/m3
nu_l = 0.294e-6;
k_l=0.6791;%thermal conductivity W/mK
sTen=0.0589; %surface tension in N/m
hlv = 2256.5e3;
% Universal constants
R_g0=8.3412; % Universal gas constant
Adis = 2e-21;

resistance = 0.5*(Tsat*(2*pi*(R_g0/M)*Tsat)^0.5)/(rho_v*hlv^2);

Rd = 0;
Rs = 0.5e-6;

% Initial conditions
del0 = (Adis/(rho_l*hlv*(Tw/Tsat - 1)))^(1/3);
per = del0/1000;
delad = del0+per;

delprimead = 0;

Pcad = Adis/delad^3;

Qad = x;   %Heat transfer per unit length

%ode solver
rspan = [Rd Rs];
y0 = [delad delprimead Pcad Qad];
odefnc = @(r,y)odefunction(r,y,sTen,Adis,hlv,rho_l,k_l,resistance,Tw,Tsat,nu_l);
options = odeset('RelTol',1e-5,'Stats','off');
[r,y] = ode15s(odefnc,rspan,y0,options);

del = y(:,1);
deldr = y(:,2);
pressure = y(:,3);
intheatflux = y(:,4);

% Curvature
K = (1/sTen)*(pressure - Adis./del.^3);

% Heat flux
heatflux = (Tw - Tsat*(1 + pressure./(rho_l*hlv)))./(resistance + del./k_l);

%Interface temperature
Tint = Tsat*(1 + pressure./(rho_l*hlv));

sfigure(2);
subplot(4,2,1)
plot(r,del)
xlabel('r')
ylabel('Film thickness \delta')
hold on
subplot(4,2,2)
plot(r,deldr);
xlabel('r')
ylabel('\delta^{\prime}')
hold on
subplot(4,2,3)
plot(r,pressure)
xlabel('r')
ylabel('Pressure \Delta P')
hold on
subplot(4,2,4)
plot(r,intheatflux)
xlabel('r')
ylabel('Heat Q [W/m]')
hold on
subplot(4,2,5)
plot(r,heatflux)
xlabel('r')
ylabel('Heat flux q')
hold on
subplot(4,2,6)
plot(r,K)
xlabel('r')
ylabel('Curvature K')
hold on
subplot(4,2,7)
plot(r,Tint)
xlabel('r')
ylabel('Local sat Temp')
hold on

pause(0.0001)

Kref = 1.2e7;
Krel = K(end)/Kref-0.1;

disp(x)
disp(Krel)
end