function [porosity_outer] = findOuterPorosityProfile(Pe,Br,zct,z,kappa,alpha)

% find critical stress
options = optimoptions('fsolve','Diagnostics','off','Display','off','MaxIterations',500);
fun = @(x)PorosityOuterEqn(x,Pe,Br,zct,z,kappa,alpha);
x0 = 0.01.*ones(size(z));
[x] = fsolve(fun,x0,options);
porosity_outer = x;

porosity_outer(z>zct) = 0;
end

