Say a universe we want to study has two finite spatial dimensions, which we will imagine to describe a square of side length two. For the diff. manifold underlying the Newtonian spacetime with which we choose to study this universe, say we equip the set $M:=\mathbb{R}^{+}_0\times[-1,1]\times[-1,1]$ with $\mathcal{O}_{\text{St.}}|_M$ and $$\mathcal{A}:=\{(M,\text{id}_M)\}$$
Our universe has, at the points corresponding to $(t,0,0)\in M$ for any $t\in\mathbb{R}^{+}_0$, a point mass of 1 unit mass. This means that on a point mass at coordinates $p=(t,x,y)$ acts a gravitational force of $$F_p:=\frac{G}{x^2+y^2}$$ Newtons. We set $G$ to $1$. The force vector on the point mass at $p$ may be described with
$$\widehat{\text{F}}_p:=(\widehat{\text{F}}^0_p,\widehat{\text{F}}^1_p,\widehat{\text{F}}^2_p)\equiv\left(0,\frac{\sqrt{1-\frac{y^2}{x^2}}}{x^2+y^2},\frac{\frac{y}{x}}{x^2+y^2}\right)$$
This determines the connection $\nabla$ (I define a connection as a map taking a vector field and a $(p,q)$-tensor field to another $(p,q)$-tensor field that satisfies the Leibnitz rule, agrees with partial differentiation for $(0,0)$-tensor fields, is $C^\infty$-linear w.r.t. the vector field input and respects addition of the tensor field input) on our manifold, as we set all the connection coefficient functions $\Gamma$ to zero except for $$\Gamma^i_{00}(p):=\widehat{\text{F}}^i_p$$
(A result of the autoparallel equation: Objects in the spacetime $M$ on which no force acts besides gravitation move along straight lines of $M$ as defined by the connection $\nabla$ if and only if we define the $\Gamma$s this way.) I want to find an injective embedding $\phi:M\rightarrow\mathbb{R}^3$ of this spacetime manifold into three-dimensional Euclidean space such that the pullback connection $\nabla^*$, coming from the Euclidean connection $\nabla^E$, is equal to $\nabla$.
The idea behind this is that we would like to see the spacetime $M$ "unstretched" - with the curvature encoded into its embedding and resulting "weird" shape - in $\mathbb{R}^3$. What is the best way to calculate what such an embedding of a spacetime would look like? And for which kind of connection carrying manifolds do such $\mathbb{R}^3$ visualizations actually exist?
My approach of finding a differential equation for an embedding $\phi$ went as follows (I would be very glad for pointing out blunders in the calculation, as otherwise solving the diff. eqs. underneath is meaningless): Per definition, we have $$\nabla^*_X\ \phi^*(s)=\phi^*\left(\nabla^E_{\phi_*(X)}\ s\right)$$
for any covectors $s\in T^1_0\mathbb{R}^3$. We will define $\phi^{-1}$ as a map from $\phi(M)$ to $M$ by $\phi^{-1}(p)=q$ iff $\phi(q)=p$. This leads to: $$\nabla^*_{e_i}\ \phi^*((\phi^{-1})^*(\in^j))=\nabla^*_{e_i}\ \in^j=\sum_q-\Gamma^{*q}_{ij}\cdot\in^q=\phi^*\left(\nabla^E_{\phi_*(e_i)}\ (\phi^{-1})^*(\in^j)\right)$$
where $e_i$ and $\in^j$ are the natural basis and induced dual basis vector fields of our chart $(M,\text{id}_M)$ respectively. We want for $\Gamma^{*q}_{ij}=\Gamma^q_{ij}$ to hold and hence get:
$$\left(\sum_q-\Gamma^q_{ij}\cdot\in^q\right)(e_k)=\left(\phi^*\left(\nabla^E_{\phi_*(e_i)}\ (\phi^{-1})^*(\in^j)\right)\right)(e_k)=\left(\nabla^E_{\phi_*(e_i)}\ (\phi^{-1})^*(\in^j)\right)(\phi_*(e_k))$$
We start to run into a notational issue here: $\left(\sum_q-\Gamma^q_{ij}\cdot\in^q\right)(e_k)$ is an element of $C^\infty(M)$, but $\left(\nabla^E_{\phi_*(e_i)}\ (\phi^{-1})^*(\in^j)\right)(\phi_*(e_k))$ is an element of $C^\infty(\phi(M))$. The above equality is therefore, technically, not correct. As the latter term stems from the definition of the covector-field-pullback $\phi^*$, the problem is that we forgot an input transformation via $(\phi(\cdot))$ at the end of it, i.e., we should write
$$\left(\sum_q-\Gamma^q_{ij}\cdot\in^q\right)(e_k)(\cdot)=\left(\nabla^E_{\phi_*(e_i)}\ (\phi^{-1})^*(\in^j)\right)(\phi_*(e_k))(\phi(\cdot))$$
We will keep this in mind and drop $(\phi(\cdot))$ again. We can calculate the left-hand side of the last result and use the Leibnitz rule for connections on the right-hand side:
$$\sum_q\delta^q_k\cdot(-\Gamma^q_{ij})=\nabla^E_{\phi_*(e_i)}\ (\phi^{-1})^*(\in^j)(\phi_*(e_k))-((\phi^{-1})^*(\in^j))\left(\nabla^E_{\phi_*(e_i)}\ \phi_*(e_k)\right)$$
This leads to:
$$-\Gamma^k_{ij}=\nabla^E_{\phi_*(e_i)}\ \in^j\bigg((\phi^{-1})^*(\phi_*(e_k))\bigg)-((\phi^{-1})^*(\in^j))\left(\nabla^E_{\phi_*(e_i)}\ \phi_*(e_k)\right)$$
$$=\nabla^E_{\phi_*(e_i)}\ \in^j(e_k)-((\phi^{-1})^*(\in^j))\left(\nabla^E_{\phi^m_{*i}e^E_m}\ \phi^n_{*k}e^E_n\right)$$
where we use the Einstein summation convention at the lower index and input of $\nabla^E$, define $e^E_i$ as the basis vector fields of $\mathbb{R}^3$ and $\phi^a_{*b}$ as the coefficient functions of the linear maps between the tangent spaces $T_pM$ and $T_p\mathbb{R}^3$ provided by $\phi_*(p)$ for each $p\in M$. Further, we get:
$$-\Gamma^k_{ij}=\nabla^E_{\phi_*(e_i)}\ \delta^j_k-((\phi^{-1})^*(\in^j))\left(\phi^m_{*i}\cdot\nabla^E_{e^E_m}\ \phi^n_{*k}e^E_n\right)$$
$$\stackrel{Leibnitz}{=}0-((\phi^{-1})^*(\in^j))\left(\phi^m_{*i}\cdot\bigg(\phi^n_{*k}\cdot\nabla^E_{e^E_m}\ e^E_n+(\nabla^E_{e^E_m}\ \phi^n_{*k})\cdot e^E_n\bigg)\right)$$
$$=-((\phi^{-1})^*(\in^j))\left(\phi^m_{*i}\cdot\bigg(\phi^n_{*k}\cdot\Gamma^{E\ a}_{mn}\cdot e^E_a+e^E_m(\phi^n_{*k})\cdot e^E_n\bigg)\right)$$
where $\Gamma^{E\ q}_{ij}$ are the connection coefficient functions of three-dimensional Euclidean space. Per definition, these functions equal zero everywhere. Hence:
$$\Gamma^k_{ij}=((\phi^{-1})^*(\in^j))\left(\phi^m_{*i}\cdot e^E_m(\phi^n_{*k})\cdot e^E_n\right)=\in^j\left((\phi^{-1})_*\bigg(\phi^m_{*i}\cdot e^E_m(\phi^n_{*k})\cdot e^E_n\bigg)\right)$$
Due to the linearity of $(\phi^{-1})_*$ we have:
$$\Gamma^k_{ij}=\in^j\left(\phi^m_{*i}\cdot e^E_m(\phi^n_{*k})\cdot (\phi^{-1})_*(e^E_n))\right)$$
By $(\phi^{-1})_{*a}^b$ we will denote the coefficient functions of the linear maps from $T_p\mathbb{R}^3$ to $T_pM$ given by $(\phi^{-1})_*$ at every point $p$ of the image $\phi(M)$. As before, instead of viewing $(\phi^{-1})_{*a}^b$ as a $C^\infty(\phi(M))$ function, as the symbol itself suggests, we will treat it as a $C^\infty(M)$ function via the transformation $(\phi(\cdot))$. Going on, we get:
$$\Gamma^k_{ij}=\in^j\left(\phi^m_{*i}\cdot e^E_m(\phi^n_{*k})\cdot (\phi^{-1})^a_{*n}\cdot e_a\right)=\delta^j_a\cdot\phi^m_{*i}\cdot e^E_m(\phi^n_{*k})\cdot (\phi^{-1})^j_{*n}$$
Which lets us conclude with:
$$\Gamma^k_{ij}=\phi^m_{*i}\cdot e^E_m(\phi^n_{*k})\cdot (\phi^{-1})^j_{*n}\ \stackrel{Einstein}{:\iff}\ \Gamma^k_{ij}=\sum_m\sum_n\phi^m_{*i}\cdot e^E_m(\phi^n_{*k})\cdot (\phi^{-1})^j_{*n}$$
The map $\phi$ we search takes values from $\mathbb{R}^3\supset(\mathbb{R}^{+}_0\times[-1,1]\times[-1,1])$ to $\phi(M)\subset\mathbb{R}^3$, which lets us exercise multivariable calculus on it:
$$\phi^a_{*b}(p)=\frac{\partial\phi^a}{\partial e_b}(p)$$ $$(\phi^{-1})^a_{*b}(\phi(p))=\frac{\partial(\phi^{-1})^a(\phi(\cdot))}{\partial e_b}(p)$$
where $\phi^a$ and $(\phi^{-1})^a$ are the $a$th component functions of $\phi$ and $\phi^{-1}$ respectively. This leads us to the following differential equations for $\phi$:
$$\left(\sum_m\sum_n\frac{\partial\phi^m}{\partial e_0}\cdot \frac{\partial}{\partial e_m}\left(\frac{\partial\phi^n}{\partial e_1}\right)\cdot \frac{\partial(\phi^{-1})^0}{\partial e_n}\right)(t,x,y)=\frac{\sqrt{1-\frac{y^2}{x^2}}}{x^2+y^2}$$ $$\left(\sum_m\sum_n\frac{\partial\phi^m}{\partial e_0}\cdot \frac{\partial}{\partial e_m}\left(\frac{\partial\phi^n}{\partial e_2}\right)\cdot \frac{\partial(\phi^{-1})^0}{\partial e_n}\right)(t,x,y)=\frac{\frac{y}{x}}{x^2+y^2}$$
and
$$\left(\sum_m\sum_n\frac{\partial\phi^m}{\partial e_i}\cdot \frac{\partial}{\partial e_m}\left(\frac{\partial\phi^n}{\partial e_k}\right)\cdot \frac{\partial(\phi^{-1})^j}{\partial e_n}\right)(t,x,y)=0$$
for any $(i,j,k)\in\{0,1,2\}^3$ where either $i\neq0$ or $j\neq0$.
Is this set of equations solvable? Can we find the solution? How to approach this? Are there conventions one should respect in their notation?
To me it seems that one big problem with this is that we don't know the codomain of our unknown $\phi$ - which is of course what I actually want to know. I am very unexperienced in diff. eqs. and have no approach to this. Is it acceptable that both $\phi$ and its inverse appear in the formula? I would also be interested to know if we can feed such complex systems of diff. equations into computational systems like Wolfram Alpha!
Edit: I noticed that, for any invertable differentiable function $f:\mathbb{R}\rightarrow\mathbb{R}$, $$f^{-1}(f(x))=x\implies\left(\frac{\partial}{\partial x}(f^{-1}\circ f)\right)(x)=1\implies\frac{\partial f^{-1}}{\partial x}(f(x))\cdot \frac{\partial f}{\partial x}(x)=1\implies\frac{\partial f^{-1}}{\partial x}(f(\cdot))=\frac{1}{\frac{\partial f}{\partial x}}(\cdot)$$
Hence, we can replace the partial deriviatives of $\phi^{-1}$ above with $$\frac{\partial(\phi^{-1})^0}{\partial e_n}(\phi(\cdot))=\frac{1}{\frac{\partial\phi^0}{\partial e_n}}(\cdot)$$
Edit 2: Unfortunately unsuccessful attempt: Assume we are really lucky and for all $(m,n)\in\{0,1,2\}^2\backslash\{(0,0)\}$ at least one of the two functions $\frac{\partial\phi^m}{\partial e_0}$, $\frac{\partial}{\partial e_m}\left(\frac{\partial\phi^n}{\partial e_2}\right)$ vanish. Then we have:
$$\frac{\frac{y}{x}}{x^2+y^2}=\left(\frac{\partial\phi^0}{\partial e_0}\cdot \frac{\partial}{\partial e_0}\left(\frac{\partial\phi^0}{\partial e_2}\right)\cdot \frac{1}{\frac{\partial\phi^0}{\partial e_0}}\right)(t,x,y)=\frac{\partial}{\partial e_0}\left(\frac{\partial\phi^0}{\partial e_2}\right)(t,x,y)$$
In more conventional notation:
$$\frac{\partial}{\partial t}\left(\frac{\partial\phi^0}{\partial y}\right)(t,x,y)=\frac{\frac{y}{x}}{x^2+y^2}$$ $$\implies\frac{\partial\phi^0}{\partial y}(t,x,y)=\int\frac{\frac{y}{x}}{x^2+y^2}\text{d}t=t\cdot\frac{\frac{y}{x}}{x^2+y^2}$$ $$\implies\phi^0(t,x,y)=\int t\cdot\frac{\frac{y}{x}}{x^2+y^2}\text{d}y$$
Freely available math software (like Wolfram Alpha) can solve this integral, and so we get:
$$\phi^0(t,x,y)=t\cdot\frac{\text{ln}(y)}{x}-t\cdot\frac{\text{ln}(x^2+y^2)}{2x}$$
If this were true, we would have one coordinate component of our embedding function $\phi$ nailed down. But it fails: If the assumption starting our current approach off holds, we have
$$\frac{\sqrt{1-\frac{y^2}{x^2}}}{x^2+y^2}=\left(\frac{\partial\phi^0}{\partial e_0}\cdot \frac{\partial}{\partial e_0}\left(\frac{\partial\phi^0}{\partial e_1}\right)\cdot \frac{1}{\frac{\partial\phi^0}{\partial e_0}}\right)(t,x,y)=\frac{\partial}{\partial e_0}\left(\frac{\partial\phi^0}{\partial e_1}\right)(t,x,y)$$
Which puts the following condition on our calculated component of $\phi$:
$$\frac{\partial}{\partial t}\left(\frac{\partial}{\partial x}\left(t\cdot\frac{\text{ln}(y)}{x}-t\cdot\frac{\text{ln}(x^2+y^2)}{2x}\right)\right)=\frac{\sqrt{1-\frac{y^2}{x^2}}}{x^2+y^2}$$
But, as simple differentiation shows, this equality is not satisfied.