(* This Mathematica program demonstrates the difference between the *)
(* non-linear response of a system and the linearized approximation *)
(* The example used is a draining spherical 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        *)

(* rad is the sphere's radius. The initial condition is h(0) = hs m.*)

(* Clear the parameter values *)

Clear[beta,rad,hs]

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

f[t_] = h'[t] + (beta*hs^(-1/2))/(Pi (2 rad - hs)) - 
       (h[t]-hs)*(beta*( rad hs^(-1/2) - 1.5 hs^(1/2) ))/
       ( Pi (2 rad hs^(1/2) - hs^(3/2))^2 )

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

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

(* Define a function which can be plotted *)

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

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

nlf[t_] = nlh'[t] + (beta*nlh[t]^(-1/2))/(Pi*(2*rad - nlh[t]))

(* Numerical solution requires parameter values *)
rad := 1
beta := 1
hs := 1/2

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

sol2 = NDSolve[{nlf[t]==0,nlh[0]==0.5},nlh[t],{t,0.0,1.3}]

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

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

(* Plot the two solutions for comparitive purposes *)

Plot[{hsolve[t],nlhsolve[t]},{t,0.0,1.25},
      AxesLabel -> {TIME,HEIGHT}]