How can I recreate the following plots from image and equations?

I'm trying to recreate these plots, so I can then manipulate them to what I want. Below is my code, the equation being modeled, and the plots I need. I can't seem to get the little kink in the profiles.

tvals = [2, 27, 3.6*365/12, 1.5*365] * 24 * 60 * 60;    %seconds
rvals = 1 : 4000;
rho=(3.6*30*24*60*60); %3.6 months
u=10^9;
u0=2*10^-5; 
v=1*10^6;
r0=2000; %Intial radius
h0=180; %Intial height
V=(3/4)*pi()*(r0^2)*h0; %volume
v0=1*10^6;
g=1.31;
tau=((3/4)^5)*(((pi()^3)*v0*(r0^8))/(g*V^3));
%% Theta 
for tidx = 1 : length(tvals)
    theta(tidx)=rho*(1-exp(-tvals(tidx)/rho));
end
%% Solving function
hz = zeros(4000, length(theta));
for ridx = 1 : length(rvals)
    r = rvals(ridx);
    for tidx = 1:length(theta)
       ta = theta(tidx);
        hz(r,tidx)=(((4*V)/(3*pi()*r0^2)).*(1./((1+(ta./tau)).^1/4)).*(1-(((r.^2)./(r0^2)).*(1./(1+((ta./tau).^1/4))))).^(1/3));
    end
end
plot(rvals, real(hz))
legend(sprintfc('%g', tvals/(24*60*60)))