Numerical method for steady-state solution to viscous Burgers' equation I am reading a paper in which a specific partial differential equation (PDE) 
on the space-time domain $[-1,1]\times[0,\infty)$ is studied. The authors are 
interested in the steady-state solution. They design a finite difference method (FDM) for the PDE. As usual, there are certain discretizations in time-space, 
$U_j^n$, that approximate the solution $u$ at the mesh points, $u(x_j,t_n)$. The 
authors conduct the FDM method on $[-1,1]\times [0,T]$, for $T$ sufficiently large 
such that
$$ \left|\frac{U_j^N-U_j^{N-1}}{\Delta t}\right|<10^{-12},\quad \forall j, $$
where $t_N=T$ is the last point in the time mesh and $\Delta t$ is the 
distance between the points in the time mesh. The approximations for the steady-state 
solution are given by $\{U_j^N\}_j$.
I wonder why the authors rely on the PDE to study the steady-state solution. As 
far as I know, the steady-state solution comes from equating the derivatives 
with respect to time to $0$ in the PDE. The remaining equation is thus an ordinary 
differential equation (ODE) in space. To approximate the steady-state solution, 
one just needs to design a FDM for this ODE, which is easier than dealing with 
the PDE for sure. Is there anything I am not understanding properly?
For completeness, I am referring to the paper Supersensitivity due to 
uncertain boundary conditions. The authors deal with the PDE $u_t+uu_x=\nu u_{xx}$, $x\in (-1,1)$, $u(
-1,t)=1+\delta$, $u(1,t)=-1$, where $\nu,\delta>0$. They employ a FDM for this 
PDE for large times until the steady-state is reached. Why not considering the 
ODE $uu'=\nu u''$, $u(-1)=1+\delta$, $u(1)=-1$, instead?
 A: You are right. Indeed, solving the ODE for steady-state (time-independent) solutions $(\tfrac12 u^2)' = \nu u''$ is rather straightforward, and we have
$$
\tfrac12 u^2 = \nu u' + k_1 \, ,
$$
where $k_1$ is an integration constant. Thus, integrating the previous ODE gives
$$
u(x) = \sqrt{2 k_1} \tanh \left( \sqrt{\frac{k_1}{2}} \frac{k_2 - x}{\nu} \right) ,
$$
where $k_2$ is an integration constant. The two integration constants are obtained by imposing the boundary conditions, and the sensitivity of the solution with respect to $\delta$ can be investigated.
With such an approach, the thing is that the time information found in the original PDE is lost. In particular, there is no more mapping with initial conditions. In numerical analysis for PDE, steady-state solutions are a good test case for the validation of time-domain numerical methods. Firstly because most of them can be computed analytically, secondly because it is quite challenging for numerical methods to converge towards the equilibrium solution as time becomes infinite. So I believe that validation was the first purpose of the authors. And then, by making the boundary data "uncertain", they remarked a sensitivity of the numerical results with respect to the boundary data.
