There are too many variables, some of them infinitely variable, to arrive at the optimal burn plan. And then you need to incorporate that burn plan into something that can be carried out with a reasonably priced computer (US $200K is "reasonably priced) than can handle high level radiation. This is not going to be a state of the art computer.
What you can do is use techniques that come close to the theoretical optimum but with constraints that keep the profile manageable by a not quite up to date computer. Mischa Kim's PhD thesis suggests using Adaptive Simulated Annealing. Markov Chain Monte Carlo is very similar to simulated annealing, and is perhaps easier to set up. Another related approach is a particle filter. All use a variant of the Metroplis algorithm, aka the Metropolis-Hastings algorithm.
Suppose you know that it will take multiple orbits from the initial orbit to the final desired orbit. In this case, you do not want the spacecraft firing continuously. You instead want the spacecraft to take advantage of the Oberth effect. In the case of very low thrust, this means firing for roughly 2/3 of the orbit, roughly centered around periapsis, plus perhaps a final finite thrust when the spacecraft nears the desired orbit.
I would do this in stages. First assume impulsive burns at periapsis, with the impulse equivalent to what could be achieved by the finite thrusters over the course of 2/3 of an orbit, with a penalty for gravity losses. Use a search technique to find a tentative multi-burn solution. This becomes the initial guess for the multiple finite burn solution. The transformation from impulsive burns to finite burns will miss the target.
So now you use similar techniques to refine the transformation. The nice thing about the Metropolis algorithm, and its variants, is that the randomness applied to tentative solutions kicks things away from local minima that in fact are lousy solutions.