A simple looking second-order DE is giving me tough time

Thank you for taking time to look at this question. Below is the equation where I need your help.

$$left ( frac{dh}{dt} right )^2 + hfrac{d^2h}{dt^2} = F_{u1} – gh $$
$F_{u1}$ and $g$ are constants.

The above second-order non-linear DE describes the following physical process:

• Consider a "T" shaped hook which pulls a lengthy heavy wire that is initially on the ground, with a constant upward force $ F_u [N]$, as shown in the below image.

• Let the height of the hook at time $t=0$ be $h=0$.

• As the wire is being pulled, due to increasing weight of the wire, the velocity of the hook would slow down, and it would finally come to rest at $t_{end}$ when static force balance is achieved.

• The variation of the hook’s height $h(t)$ is desired to be known.

To determine $h(t)$, Newton’s second law is used as follows:
$$ F_{net} = frac{d (mv)}{dt}, text{where} v=frac{dh}{dt} $$
$$ F_u – F_{wire} = frac{d (rho A h frac{dh}{dt})}{dt}$$
$$ F_u – rho A h g = rho A frac{d (h frac{dh}{dt})}{dt}$$
$$ frac{F_u}{rho A } – h g = left ( frac{dh}{dt} right )^2 + hfrac{d^2h}{dt^2} $$
$$left ( frac{dh}{dt} right )^2 + hfrac{d^2h}{dt^2} = F_{u1} – gh $$
$$ text{where } F_{u1}=frac{F_u}{rho A}$$

As closed-form analytical solution is not available, finite difference is used as follows:

I tried to solve it using the below code. As the above equation is second-degree polynomial in $h^{t+Delta t}$, there will be two solutions. The solutions will be of opposite sign until the equilibrium height is reached. So the code selects only the positive value of the root, which has physical meaning.

I took the initial height of the wire, h(2), as 0 m. h(1) is treated as a ghost point in time, i.e. Time t = 0 at h(2), t = dt at h(3), t = 2*dt at h(4) and so on. Here, h(1) is at t = – dt, and hence h(1) = h(2). I get the solution as shown below, where the red line indicates the height from static force balance.

As you could see from the solution plot, the numerical solution reaches an equilibrium height that is considerably higher than the height from static force balance. I tried reducing the time step, but the difference remains the same.

Could you please shed some light ? Thank you.

clc;
clear all;


rho     =1000; 
radius  =3; 
area    =pi*(radius/1000)^2; 
g = 9.81; 


F_u=1; 

F_u1 = F_u/(rho*area); 

%static force balance
h_static = F_u/(rho*area*g);

dt = 0.00001; 

h(1) = 0;
h(2) = 0; 
%Model equation Ax^2 + Bx + C = 0

dt_steps = 8/dt;
for i =2:1:dt_steps
A=1;
B=h(i);
C = (h(i))^2-(h(i)*h(i-1))+((F_u1-(g*h(i)))*dt^2);
polynomial_roots=roots([A -B -C]);

if (length(polynomial_roots(polynomial_roots>0))==1)
h(i+1)=polynomial_roots(polynomial_roots>0);
else 
    break
end
end

scatter(((1:length(h))-2)*dt,h);
hold on;
plot([0 (length(h)-2)*dt],[h_static h_static],'-r');
xlim([0 (length(h)-2)*dt]);
xlabel('Time [s]')
ylabel('Height [m]')
title(['dt = ' num2str(dt),' s'])