(* This Mathematica program demonstrates the difference between the *)
(* non-linear response of a system and the linearized approximation *)
(* The example used is a draining conical tank in which the         *)
(* outflow is proportional to the square root of the liquid level   *)

(* The parameters used are Fi [=] m^3/hr, beta [=] m^2.5/min        *)

(* Theta is the cone's angle. The initial condition is h(0) = hs m. *)

(* Clear the parameter values *)

Clear[beta,theta,hs]

(* Define a function which is the linear o.d.e. *)

f[t_] = h'[t] + (beta*hs^(-3/2))/(Pi (Tan[theta])^2) - 
                (h[t]-hs)*(3*beta*hs^(-5/2))/(2*Pi*(Tan[theta])^2)

(* Solve the linear o.d.e. with the initial tank height *)
(* set to hs m. *)

 sol2 = DSolve[{f[t]==0,h[0]==hs},h[t],t]

(* Define a function which can be plotted *)

hsolve[t_] = h[t] /. sol2[[1]]

(* Define a function which is the non-linear o.d.e. *)

nlf[t_] = nlh'[t] + ( beta/(Pi (Tan[theta])^2) )*nlh[t]^(-3/2)

(* Numerical solution requires parameter values *)
Fi := 0
theta := 0.52
beta := 1
hs := 5

(* Solve the non-linear o.d.e. with the initial tank *)
(* height set to hs m. *)

sol = NDSolve[{nlf[t]==0,nlh[0]==hs},nlh[t],{t,0,23}]

(* Convert the numerical solution to a function which *)
(* we can plot *)

nlhsolve[t_] =  nlh[t] /. sol[[1]]

(* Plot the two solutions for comparitive purposes *)

Plot[{hsolve[t],nlhsolve[t]},{t,0,23}]