Can one introduce magnetic monopoles without Dirac strings?
If you have a source of radial magnetic field $B\sim Q_M/r^2$, then one may prove that the vector potential $\vec A$ can't be single-valued. It's because $\vec B={\rm curl}\vec A$ for a well-defined $\vec A$ automatically satisfies ${\rm div}~\vec B=0$. However, $Q_M/r^2$ has a curl proportional to the delta-function at the origin.
Still, this delta-function vanishes everywhere except for a Dirac string (and the space minus the Dirac semiinfinite string is simply connected), so with the Dirac string, $\vec A$ may be defined everywhere. $\vec A$ still changes under the loop around the Dirac string. In this way, the magnetic monopole is replaced by a very long magnetic dipole. Two poles are connected with a thin solenoid and one of the monopoles is sent to infinity and becomes irrelevant. The Dirac string, i.e. a very thin solenoid, becomes unobservable as well, even for interference experiments (as long as the confined magnetic flux is properly quantized).
The arguments above are waterproof and you can't circumvent them. So if you're asking whether there is a way to introduce a magnetic monopole so that the vector potential would be single-valued, the answer is a resounding No, much like if you ask if it is possible to introduce the number 4 so that it isn't equal to 2+2.
However, one may try to search for solutions to similar, less singular problems. In theories with Higgses, one may "dilute" the delta-function a little bit and find non-singular solutions of Yang-Mills theories with Higgs fields, the so-called
http://en.wikipedia.org/wiki/%27t_Hooft%E2%80%93Polyakov_monopole
't Hooft-Polyakov monopole that is non-singular but indistinguishable from the Dirac monopole when you're very far from the center of the solution, relatively to its characteristic length scale. This solution has various generalizations, too.
There is also the geometric approach to a Dirac-type magnetic monopoles with a singular core, pioneered by Wu and Yang, which avoids the Dirac string by introducing non-trivial topology, fiber bundles, and locally defined gauge potentials. It leads to the same physical predictions, e.g. charge quantization, as Dirac's method.
(The above should not be conflated with a 't Hooft-Polyakov-type magnetic monopoles with a regular core, cf. Lubos Motl's answer and this Phys.SE post.)
References:
T.T. Wu and C.N. Yang, Concept of non-integrable phase factors and global formulation of gauge fields, Phys. Rev. D 12 (1975) 3845.
T.T. Wu and C.N. Yang, Dirac Monopole without Strings: Classical Lagrangian theory, Phys. Rev. D 14 (1976) 437.
M. Nakahara, Geometry, Topology and Physics, 1990.
Introducing magnetic charge into Maxwell equations is not a problem at all, and it does not require any strings etc. Moreover, it makes Maxwell equations symmetric w.r.t. magnetic and electric fields/charges. The equations are as follows:
\begin{split} \mathop{\mathrm{curl}} E + \frac{\partial H}{\partial t} &= -J_m \\ \mathop{\mathrm{curl}} H - \frac{\partial E}{\partial t} &= J_e \\ \mathop{\mathrm{div}} E &= J_e^0 \\ \mathop{\mathrm{div}} H &= J_m^0 \end{split}
However, introduction of magnetic charge leads to non-zero divergence of magnetic fields, hence making it impossible to represent the magnetic field as a curl of vector potential. As you stated in your question, Dirac introduced a singularity to preserve gauge potentials description of the phenomenon.
Use of gauge potentials is justified by an excellent fit of QED predictions and experimental data in all aspects that are not related to identification of particle's mass and charge values from the theory. The straightforward calculation leads to infinities, and re-normalization procedure does not help with identification of mass and charge values.
Elementary particles that were observed experimentally are known to have zero magnetic charge but non-zero magnetic moment. Hence there might be a non-zero distribution of magnetic charge density "inside" the particles, if you admit that particles are not point-like. This approach can be developed using Gordon decomposition of the vector current constructed from Dirac spinors. See for instance here.