In its broadest form, Trajectory Optimization is simply finding the best path (and control values to produce this path) that satisfy the dynamics. The dynamics, in this case, are either the orbital dynamics (e.g. kepler's 2 body problem with proper ephemeris), rocket ascent dynamics or ballistic/lifting reentry dynamics. The cost function can be created in many different ways, depending on what you want in the end. For example, you may want to optimize the dynamics to produce a particular state at the end of a fixed time interval. Or you may want to minimize the fuel used in the trip (which maximizes the payload mass). Or you may want to minimize the time it takes to reach a particular end state using the least amount of (throttleable/vectored thrust).
Once you establish what your dynamic equations are and what your "cost function" are, you can translate that into a set of equations which you can solve. This process is called "transcription" in trajectory optimization nomenclature. Broadly speaking, there are two ways to transcribe the problem. The "indirect" method suggests that you should formulate the optimality conditions first analytically, then discretize those conditions over some number of timesteps N. This will result in a non-linear system of equations with $M\cdot N$ variables, where M is the number of variables to describe a state at a given time (e.g. the number of differential equations); plus $P\cdot N$ control variables, where P is the number of controls that are applied in each timestep. For example, in an 2-body problem in cartestian coordinates, you will have four state variables at each timestep (x-position, x-velocity, y-position, and y-velocity). The control may be a throttleable thrust which always points in the velocity direction. So, there would be only 1 control at each timestep. Often times, this can be tedious to do since a lot of the work has to be done by hand.
The alternative is referred to as the "direct" method. In this case, you would discretize the equations first and then formulate the conditions for optimality. So instead of formulating the conditions on an infinite space (indirect method), you're formulating them only for the finite dimensional space (a finite number of state + control variables). The advantage of this method is that the problem now becomes a non-linear constrained optimization problem (which is easier to program directly and solve using software packages).
Under the umbrella of direct methods, there are couple of different ways to discretize the problem. One way is called collocation, where we force the dynamics to agree at certain points in time (while assuming some sort of interpolation scheme between these points, such as linear or quadratic interpolation). The other approach is called shooting, which is kind of like a trial and error approach where the discrepancy (a.k.a. the defect constraint) at each timestep is discretized and also minimized. I haven't used shooting methods for trajectory optimization, but I have used collocation method so I will address it directly.
The key trick to the collocation method is knowing how to compute the non-linear equality constraints using the dynamics of the problem. The key idea is that on each time step, one must constrain the predicted value from the dynamics to be equal to the true value at the end of the timestep. Suppose we are considering the ith timestep, which has state variables $x_i$ and $x_{i+1}$ respectively at timesteps $t_i$ and $t_{i+1}$. Let's also suppose that the dynamic equation can be written in the form $\dot{x}=f$. If we integrate both sides of this equation with respect to time over the interval from $t_i$ to $t_{i+1}$, we get $$x_{i+1}-{x_i}=\int_{t_i}^{t_{i+1}}f dt.$$
One can discretize the integral in many ways. For simplicity, let's choose something like the trapezoid rule. Then we can obtain $$x_{i+1}={x_i}+\frac{f_{i}+f_{i+1}}{2} \Delta t.$$
Let's define the value predicted by the dynamics as $x_{i+1}^p={x_i}+\frac{f_{i}+f_{i+1}}{2} \Delta t.$ The actual value at the end of this timestep is $x_{i+1}$. So, to impose a dynamic constraint, we simply set the predicted value equal to the true value at the end of this time interval. Mathematically, we write: $x_{i+1}^p-x_{i+1} =0$. We write it this way because most nonlinear programming packages require you to write equality constraints with the RHS equal to zero.
That's the majority of the hard work. You simply need to create a function which computes a vector containing all of these constraints at each timestep simultaneously. Similarly, you will have to write your objective function in terms of the discrete variables $x_i$, which may include a discretized integral over time and/or a function of the end state.
Once you have the objective function and the equality constraint function (both of which are functions of the state variables and controls at all timesteps), you will also have to provide lower and upper bounds on each of the state and control variables at all timesteps. These bounds are usually physically motivated, like mass must be positive or the orbit radius cannot exceed 2 astronomical units, etc... Additionally, you often have to provide an initial guess that satisfies the lower and upper bounds. Sometimes, a trajectory optimization problem may be sensitive to your initial guess so it is a good idea to choose it as close to what you think the true solution should be as possible (e.g. if you have any physical intuition into the problem, now is the time to use it).
Once you have all this information, you can just feed it to a nonlinear programming solver (e.g. Matlab's fmincon, ipopt, etc...). There are usually different methods that you can use to solve nonlinear constrained optimization problems, such as sequential quadratic programming (SQP), active set methods (ASM), or interior point methods (IP). Of the methods I've used in trajectory optimization for orbital mechanics, I have found interior point methods to be the fastest and most robust. Regardless of which non-linear programming method you choose, make sure that when you check for convergence, check both feasibility (how much constraint violation are we incurring) and first-order optimality (how close are we to being at the bottom of a bowl/valley). Both values need to be small in order to be considered a "converged" solution.
Before you start programming a trajectory optimization code for orbital mechanics, I highly recommend trying a simpler problem with a known solution. The so-called bang-bang control problem is a good example of a 1D optimal control problem involving 1D motion. IN this problem, one seeks the optimal acceleration required to reach a certain distance in the minimum time such that your final velocity is also zero. This is a classic problem with lots of physical intuition (and even an analytical solution). Intuitively, the best control would be to accelerate as fast as you can toward your goal destination up until a certain point and then decelerate as fast as possible until you just barely reach your destination point with a velocity of zero.
You can see how this bang-bang problem can be transcribed using the direct collocation method (complete with matlab code) on Sam Pfrommer's blog:
http://sam.pfrommer.us/tutorial-direct-collocation-trajectory-optimization-with-matlab
For additional information on trajectory optimization in general, with more complicated examples such as the cart-pole control problem, I suggest reading/watching Matthew Kelley's tutorial on this topic:
https://epubs.siam.org/doi/pdf/10.1137/16M1062569
