This is a problem in modeling in hydraulic fracturing field. It's quite long so hopefully someone can patiently read and help me. 
The equation numbers are match those of the reference paper by Adachi, Detournay and Peirce.
Mathematical model:
The fracture profile $w(z)$ can be calculated explicitly for a given height $h$:
$$ \begin{split} w &=\frac{2}{E'}\left(p-\Delta \sigma \sqrt{h^2-4z^2}\right) + \frac{4 \Delta \sigma}{\pi E'} \sqrt{h^2-4z^2} \arcsin \frac{H}{h} \\ & \; - z\ln\left|\frac{H\sqrt{h^2-4z^2}+2z\sqrt{h^2-H^2}}{H\sqrt{h^2-4z^2}-2z\sqrt{h^2-H^2}}\right|+\frac{H}{2}\ln\left|\frac{\sqrt{h^2-4z^2}+\sqrt{h^2-H^2}}{\sqrt{h^2-4z^2}-\sqrt{h^2-H^2}}\right| \end{split}\label{7}\tag{7}$$ where $$ 2\arcsin\frac{H}{h}-\left(1-\frac{p}{\Delta\sigma}\right)\pi=\sqrt{\frac{2\pi}{H}}\frac{K_{IC}}{\Delta\sigma} \sqrt{\frac{H}{h}}.\label{8}\tag{8} $$ From \eqref{7} & \eqref{8} it follows that $$ \overline{w}=\overline{w}(h),\quad p=p(h),\label{9}\tag{9} $$ where: $$ \overline w=\frac{1}{H} \int\limits_{-h/2}^{h/2} wdz,\quad \overline q=\frac{1}{H} \int\limits_{-h/2}^{h/2} qdz\label{1}\tag{1} $$ Poiseuille's law $q(z,t)$: $$ q=-\frac{w^3}{12 \mu} \frac{\partial w^3}{\partial x},\qquad-h/2<z<h/2 \label{10}\tag{10} $$ After integration over $h$ yields: $$ \overline q=-\frac{\psi' w^3}{\mu} \frac{\partial w^3}{\partial x}\label{11}\tag{11} $$ where $\psi'$: shape factor
Conservation law:
$$ \frac{\partial \overline q}{\partial x}+ \frac{\partial \overline w}{\partial t}+\frac{C'}{\sqrt{t-t_0(x)}}=0\label{12}\tag{12} $$
where $C'=2C_l$, with $C_l\equiv$ leak-off coefficient
Boundary conditions:
At $x=l(t)$:
$$
\overline{w}(l,t)=\overline{q}(l,t)=0,\quad h(l,t)=H.\label{13}\tag{13}
$$
At the inlet:
$$
\overline{q}(0,t)=\frac{Q_0}{2H}.\label{14}\tag{14}
$$
Scaling:
$$
\begin{matrix}
\xi =x/l, & \zeta=z/h_*,
\\ \tau=t/t_* , & \gamma=l/l_* ,
\\ \lambda=h/h_*, & \Pi=p/p_* ,
\\ \Psi=q/q_* ,& \overline \Psi=\overline q/q_* ,
\\ \Omega=w/w_* ,& \Omega =\overline w/w_*
\end{matrix}\label{15}\tag{15}
$$
and
$$
\begin{matrix}
l_*=\dfrac{\pi H^4 \Delta \sigma^4}{4E'^3 \mu Q_0}, &
h_*=H,\\
t_*=\dfrac{\pi^2 H^6 \Delta \sigma^5}{4 \mu E'^4 Q_0^2}, &
p _*=\Delta \sigma, \\
q=\dfrac{Q_0}{2H}, &
w=\dfrac{\pi H \Delta \sigma}{2E'}=M_0 \Delta \sigma
\end{matrix}\label{16}\tag{16}
$$
Scaled governing system:
$$
\begin{eqnarray}
\frac{\partial \overline \Omega}{\partial \xi}=-\frac{\gamma \overline \Psi}{\Upsilon(\overline \Omega,\kappa) \overline {\Omega}^3}=F_1(\xi,\overline \Omega,\overline \Psi,\gamma)\label{55}\tag{55}\\
\frac{\partial \overline \Psi}{\partial \xi}=-\frac{\xi \gamma \dot{\gamma} \overline \Psi}{\Upsilon(\overline \Omega,\kappa) \overline {\Omega}^3}-\gamma \dot{\overline \Omega}-\frac{C}{\sqrt{\tau-\tau_0(\gamma \xi)}}=F_2(\xi,\overline \Omega,\overline \Psi,\gamma)\label{56}\tag{56}\\
\tau=\gamma \int_0^1 {\overline{\Omega}(\xi \tau)}d\xi+2C \gamma \int_0^1 {\sqrt{\tau-\tau_0(\gamma \xi)}}d\xi \label{57}\tag{57}
\end{eqnarray}
$$
Boundary conditions at the inlet and at the leading edge:
$$
\begin{cases}
\overline{\Psi}(0,\tau)=1,\\
\overline{\Omega}(1,\tau)= \overline{\Psi}(1,\tau)=0
\end{cases}\label{39}\tag{39}
$$
The nonlinear governing system was derived using backward finite difference, composite Simpson law, Hermite... and then solved using Newton Jacobian (matlab codes is provided here
Question
Now I want the inlet injection rate $\overline{q}(0,t)$ changes in time so I can run simulation using the provided codes in the case of a variable rate. Originally, $\overline{q}(0,t)$ is constant so the scaled injection rate is $\overline{\Psi}(0,\tau)=1$. If I let it changes in time, what should be the form of the injection rate at the inlet $\overline{q}(0,t)$ and scaled injection rate $\overline{\Psi}(0,\tau)$ accordingly?