"Is there a change I can make to the code which will get it working?"
No, there is not. Your substitution of -exp(-y) for theta in only one place messes up the equation, including the boundary condition.
The original equation is
diff(theta(zeta), zeta, zeta) = -2*(diff(theta(zeta), zeta))/zeta - theta(zeta)
and your revised version is
diff(y(x), x, x) = -2*(diff(y(x),x))/x + exp(-y(x))
Notice that the -theta(zeta) has changed to +exp(-y(x)) which is a change of sign as well as a change of variables.
You can make the substitution into the theta/zeta equation of
theta(zeta) = -exp(-Y(zeta))
followed by zeta = x
This would get you
(diff(Y(x), x, x))*exp(-Y(x))-(diff(Y(x), x))^2*exp(-Y(x)) = -2*(diff(Y(x), x))*exp(-Y(x))/x + exp(-Y(x))
which has the proper final term but now has the extra exp(-Y(x)) terms. We can divide through by that extra term to hope to recover the front terms:
-(diff(Y(x), x))^2+diff(Y(x), x, x) = (-2*(diff(Y(x), x))+x)/x
but instead we get different front terms and the +exp(-Y(x)) has become +1
If you had stuck to the original form of the equation with the obvious substitutions,
diff(y(x), x, x) = -2*(diff(y(x),x))/x - y(x)
then dsolve() can handle that by itself in a general form, or can give you a useful answer if you add the initial condition for y(x)=0, but it does seem to have difficulty if you add the other initial condition.
