Dynamical systems occur in many areas of science, especially fluid dynamics. One is often interested in examining the structural changes in dynamical systems, and these are often related to the appearance or disappearance of solution trajectories connecting one or more stationary points. Homoclinic orbits are trajectories that connect a hyperbolic saddle-type stationary point to itself, and often arise as the limiting case of periodic solutions. In global bifurcation analysis, the occurrence of homoclinic orbits is closely tied to chaotic behaviour in dynamical systems. Numerical computation of homoclinic orbits requires solving a boundary value problem (BVP) over a doubly infinite interval. There are three ways that this can be done. The first is to truncate the problem to a large finite interval and create appropriate projected boundary conditions to ensure that the solution over the infinite interval is approximated well. The second is to transform the problem to one over a finite interval via an exponential change of independent variable. The third is to use a method appropriate for an infinite interval by discretising the problem using a set of approximating functions appropriate to an infinite interval. In this thesis we will compare algorithms of all three types, and test our implementations of them on several test problems. • The first algorithm we will examine will be the original time method of Beyn. This method truncates the infinite interval to a finite one, and produces projected boundary conditions that force the endpoints of the trajectory to lie in the appropriate invariant subspaces of the stationary point. This algorithm also includes a method for determining whether or not the truncated interval is large enough to obtain a solution accurate to a required tolerance. • The second algorithm is called time + subspace, and is a variation on the original time method that fixes the distance between the stationary point and the endpoint of the solution, but allows the interval to vary – the reverse of which was true for original time. • The third algorithm is of the second type and is based upon the arclength parameterisation of the orbits. This method uses the arclength of the trajectory as the independent variable to transform the problem to a finite interval. However, the BVP formulated in this way can have a singularity at the end of the domain, and thus a special collocation method is required to handle this. • As this method does not perform well on orbits exhibiting ˇ Silnikov behaviour, a variation of this method called the partial arclength method which uses a different exponential transformation of the independent variable near the stationary point is presented to address the deficiency. • The final algorithm presented in this thesis uses Laguerre polynomials to compute points on the stable and unstable manifolds of the stationary point by solving a BVP over a semi-infinite interval. These are then used to produce boundary conditions for the truncated BVP that ensure that the endpoints are on the appropriate manifold and not just the appropriate subspace.