prediction of alpha given maximum height.

The code which I posted below calculates the maximum height.This height depends on the entrainmnet constant(alpha) which I used in the last function. Now my motive is, let's say the max height is given and I want to calculate alpha. I have tried different ways none of them not working well. Thanks for the help.

  zspan=[0,400];
v0mat = [1 0.01 1];
zsol = {};
v1sol = {};
v2sol = {};
v3sol = {};
for k=1:size(v0mat,1)
v0=v0mat(k,:);
[z,v]=ode45(@rhs,zspan,v0);
[~, Nsqr]=density1(z);
zsol{k}=z;
v1sol{k}=v(:,1);
v2sol{k}=v(:,2);
v3sol{k}=v(:,3);
end
for r=1:length(ZH)
q(r)=r;
end  
for k2 = 1:length(ZH)
Mzsol(k2)= max(ZH{k2});
end
%Plot of maximum height. 
figure()
scatter(q,Mzsol,'g')
xlabel('q')
ylabel('Maximum height')
function [rho, Nsqr]=density1(z)
%%Data for desnity with respect to depth 
zval= [2 3 5 7 10 15 20 25 30 40 50 60 70 80 90 100 125 150 160 175 200 225 250 275 300 325 350 375 400];
rho1 = [17.2731684 17.1649375 21.43455647 22.4140625 23.86332207 24.3746967 24.70487685 24.6003125 24.8933125 25.42772826 26.03220776 26.439625 26.8151875 26.86830797 27.1949375 27.34406944 27.5551875 27.728625 27.23423729 27.88542857 27.752249049 28.1025 28.2415 28.37 28.05366667 28.6565 28.7755 28.898 29.013]; 
rho0=29;g=9.8;
zvala=400-zval;
zvalue=fliplr(zvala);
rho2=fliplr(rho1);
rho3=smoothdata(rho2,'lowess',5);
rho=interp1(zvalue,rho3,z);
rho4=interp1(zvalue,rho3,z+0.1);
N=(-g./rho0).*(rho4-rho);
Nsqr= sqrt(N);
end
function parameters=rhs(z,v)
[~, Nsqr]=density1(z);      
%Entrainmnet constant
alpha=0.116; 
        db= 2*alpha-(v(1).*v(3))./(2*v(2).^2);
        dw= (v(3)./v(2))-(2*alpha*v(2)./v(1));
        dgmark= -Nsqr.*Nsqr-(2*alpha*v(3)./v(1)); 
parameters=[db;dw;dgmark];
end

[Dereje]

I am not sure if shooting method helps!