A model, composed of the nonlinear Cahn-Hilliard and Flory-Huggins theories, is used to numerically simulate the phase separation and pattern formation phenomena of a representative oligomer solution when it is quenched into the unstable region of its binary phase diagram. This model incorporates initial thermal concentration fluctuations, which are computed using a Monte Carlo technique, and zero mass flux and natural nonperiodic boundary conditions. A computational method is presented for the numerical solution of this previously unsolved realistic model. This method uses Galerkin finite elements with Hermitian bicubic interpolants to discretize space and a finite difference method to discretize time. In addition, it uses a first-order implicit Euler predictor-corrector method for time integration, and an adaptive time-step controller. The numerical results replicate experimental observations when binary solutions are quenched into the unstable region; i.e., the interconnected structure forms for critical quenches, and a droplet-type morphology forms for off-critical quenches. A dimensionless diffusion coefficient D is identified, which controls the amount and rate of phase separation; as D increases, the amount and rate of phase separation increase as well. Lastly, the initial concentration affects the lag time for macroscopic phase separation to first occur. For a given temperature that is less than the critical temperature, this time is smallest at the critical concentration and increases to infinity at the two spinodal concentrations.