Solving systems of partial differential equations

I have been trying to find a way for a while now to simulate a few different systems of PDEs. What I need to be able to do is change both the initial conditions and the system of equations, without having to go through and make changes to how I am solving the system every time I make a change to what the system is that I am solving.

The systems of equations is closely related to the maxwell equations and this suggested the idea of using FDTD, I have since opted against this idea for the simple reason that I am not sure how to set up the boundary conditions for what I am doing and I am concerned that when I change the equations that I am simulating the needed boundary conditions would change in ways that would be difficult to solve for and lead to problems for the simulation that I can’t actually know the effect of without solving the equations.

This led me to the idea of using NeuralPDE, in fact I got as far as writing a simple case or two when I got kind of concerned about things like guaranteeing convergence. In particular I am concerned that I will get solutions that are valid only in special cases as an example, where one of the functions has been set to a constant.

At this point I looked around to see if there was an obvious choice for doing something like this mathematically, what I have found is the Homotopy Analysis Method. It looks like this is just what I have been looking for from a mathematical stand point. Unfortunately I have not been able to find a Homotopy solver for PDEs written in Julia. Is there such a thing and I am just missing it, or is there another choice out there that I should consider, or am I better off if I write one?