Stochastic differential equation (SDE) is a differential equation in which some of the terms and its solution are stochastic processes. Stochastic differential equations can be formulated from given ordinary differential equations by introducing stochastic perturbations in it. SDEs play a central role in modeling physical systems like finance, biology, engineering to mention some. In modeling process, the computation of the trajectories (sample paths) of solutions to SDEs is very important. However, the exact solution to an SDE is generally difficult to obtain due to non-differentiability character of realizations of the Brownian motion. There exist approximation methods of solutions of SDEs. The solutions will be continuous stochastic processes that represent diffusive dynamics, a common modeling assumption for financial, biological, physical and environmental systems. This thesis presents a comprehensive survey on numerical methods for stochastic differential equations. Standard numerical methods, Euler Maruyama, stochastic Euler, stochastic Runge Kutta methods are well described. These methods have certain limitations and do not behave well in certain scenarios. A reliable stochastic non-standard finite difference (SNSFD) scheme is proposed which remains consistent for all choice of parameters. In this thesis, we described formulation of nonparametric perturbation techniques for stochastic models. The formulation is based on epidemiological model of humans, animals and plants. Also, this idea is extended for computer virus transmission and smoking models. A comparison of proposed stochastic non-standard finite difference (SNSFD) method with standard stochastic numerical methods is also presented. Explicit methods, such as Euler Maruyama, stochastic Euler, stochastic Runge-Kutta methods are widely used in solving systems of stochastic differential equations. It is well-known that solving stochastic differential equations with explicit finite-difference schemes such as Euler Maruyama, stochastic Euler, stochastic Runge-Kutta methods can result in contrived chaos and non-physical oscillations caused by numerical instability for certain values of the discretization parameters. Such scheme-dependent numerical instabilities can be avoided by using small time steps but the additional computing cost incurred when examining the long-term behavior of a dynamical system. To avoid unnatural chaos and other scheme dependent numerical instabilities, implicitly-driven explicit scheme with certain additional desired properties are generally preferred. The numerical integration of systems of stochastic differential equations over very long-time intervals requires the use of time steps which are the largest possible, keeping in mind accuracy and stability. In the present research work, we developed and investigated such reliable numerical method for the solution of stochastic models, which remains consistent with the continuous dynamical systems. This method will numerically analyze the behavior of solution of the models, stability analysis of the steady states and threshold criteria for the physical systems. The proposed method could be used with arbitrarily large time steps, thus making them more economical to use when integrating over long time periods and could be restore all the essential properties like dynamical consistency, positivity and boundedness of the corresponding dynamical systems.