Since you already have an expression for $x'$, namely, $x'=f(x,p)$, we can simply take the differential of that expression:

$$\delta x'=\frac{\partial f}{\partial x}\delta x+\frac{\partial f}{\partial p}\delta p$$

It's not entirely clear what your procedure actually is, but I'll write this answer assuming that it's something like the following (let's assume for simplicity of notation that there are three spatial variables $x,y,z$ and three parameters $a,b,c$):

You have a system of differential equations:

$$x'=f(x,y,z,a,b,c)$$ $$y'=g(x,y,z,a,b,c)$$ $$z'=h(x,y,z,a,b,c)$$

This system of differential equations has a solution:

$$x=F(t,a,b,c)$$ $$y=G(t,a,b,c)$$ $$z=H(t,a,b,c)$$

which depends both on time and on the parameters $a,b,c$.

You also have a numerical solver that can, for a particular value of $a,b,c$, generate $F,G,H$ from $f,g,h$.

To determine the best fit, you compare $F(t)$, $G(t)$ and $H(t)$ to some input data $x_0(t)$, $y_0(t)$, and $z_0(t)$, varying $a,b,c$ to try to minimize some measure of the difference between the numerical solution and the input data. From there, you get your optimal parameters.

So what you're actually trying to do is find the uncertainty in a function of three variables and three functions of four variables: namely, you're trying to find the uncertainty in the following quantity:

$$f(x(t,a,b,c),y(t,a,b,c),z(t,a,b,c),a,b,c)$$

which has both explicit and implicit dependence on the parameters.

Generally, we cannot assume that the parameters are uncorrelated. If you're using any kind of decent fitting software, you should be able to get a covariance matrix $C$ which contains information about the correlations, defined as follows:

$$C=\begin{pmatrix}\sigma_a^2&\sigma_{ab}&\sigma_{ac}\\\sigma_{ba}&\sigma_b^2&\sigma_{cb}\\\sigma_{ca}&\sigma_{cb}&\sigma_c^2\end{pmatrix}$$

where, for example, $\sigma_a^2$ is the variance of $a$, and $\sigma_{ab}$ is the covariance between $a$ and $b$. From there, we can get the uncertainties $\sigma_x$, $\sigma_y$, and $\sigma_z$ using the solution $F,G,H$, the elements of this matrix, and the timestep $\sigma_t$:

$$\sigma_x^2=\left(\frac{\partial F}{\partial t}\right)^2\sigma_t^2+\left(\frac{\partial F}{\partial a}\right)^2\sigma_a^2+\left(\frac{\partial F}{\partial b}\right)^2\sigma_b^2+\left(\frac{\partial F}{\partial c}\right)^2\sigma_c^2+2\frac{\partial F}{\partial a}\frac{\partial F}{\partial b}\sigma_{ab}+2\frac{\partial F}{\partial a}\frac{\partial F}{\partial c}\sigma_{ac}+2\frac{\partial F}{\partial b}\frac{\partial F}{\partial c}\sigma_{bc}$$

and similarly for $y$ and $z$, substituting $G$ and $H$ respectively for $F$.

The part I'm less sure about comes in here: we also cannot assume that the uncertainties in $x,y,z$ are uncorrelated, so we need to propagate the covariances as well. I believe this is done in the following way:

$$\sigma_{xy}=\left\vert\frac{\partial F}{\partial t}\right\vert\left\vert\frac{\partial G}{\partial t}\right\vert\sigma_t^2+\left\vert\frac{\partial F}{\partial a}\right\vert\left\vert\frac{\partial G}{\partial a}\right\vert\sigma_a^2+\left\vert\frac{\partial F}{\partial b}\right\vert\left\vert\frac{\partial G}{\partial b}\right\vert\sigma_b^2+\left\vert\frac{\partial F}{\partial c}\right\vert\left\vert\frac{\partial G}{\partial c}\right\vert\sigma_c^2+\left(\frac{\partial F}{\partial a}\frac{\partial G}{\partial b}+\frac{\partial G}{\partial a}\frac{\partial F}{\partial b}\right)\sigma_{ab}+\left(\frac{\partial F}{\partial a}\frac{\partial G}{\partial c}+\frac{\partial G}{\partial a}\frac{\partial F}{\partial c}\right)\sigma_{ac}+\left(\frac{\partial F}{\partial b}\frac{\partial G}{\partial c}+\frac{\partial G}{\partial b}\frac{\partial F}{\partial c}\right)\sigma_{bc}$$

and, again, similarly for the other pairs of variables.

So now we have our variances and covariances for all of our variables in the main function, which means we can finally evaluate $\sigma_{x'}$:

$$\sigma_{x'}^2=\left(\frac{\partial f}{\partial x}\right)^2\sigma_x^2+\left(\frac{\partial f}{\partial y}\right)^2\sigma_y^2+\left(\frac{\partial f}{\partial z}\right)^2\sigma_z^2+2\frac{\partial f}{\partial x}\frac{\partial f}{\partial y}\sigma_{xy}+2\frac{\partial f}{\partial x}\frac{\partial f}{\partial z}\sigma_{xz}+2\frac{\partial f}{\partial y}\frac{\partial f}{\partial z}\sigma_{yz}+\left(\frac{\partial f}{\partial a}\right)^2\sigma_a^2+\left(\frac{\partial f}{\partial b}\right)^2\sigma_b^2+\left(\frac{\partial f}{\partial c}\right)^2\sigma_c^2+2\frac{\partial f}{\partial a}\frac{\partial f}{\partial b}\sigma_{ab}+2\frac{\partial f}{\partial a}\frac{\partial f}{\partial c}\sigma_{ac}+2\frac{\partial f}{\partial b}\frac{\partial f}{\partial c}\sigma_{bc}$$

and similarly for $\sigma_{y'}$ and $\sigma_{z'}$, substituting $g$ and $h$ for $f$. It should be fairly straightforward to see how to generalize this to higher numbers of parameters and higher dimensions of solution vector.