3D parametric equations for an elliptical orbit, using inclination angle

Question

I am modelling planetary orbits using parametric equations, with the Sun at (0,0).

In 2D, it's quite straightforward ($a$ is the semi-major axis length and $e$ is eccentricity)

$$ x\left(t\right)=a\ cos\left(t\right)+ae,\\ y\left(t\right)=a\sqrt{1-e^2}\ sin\left(t\right)\\$$

However, I'm struggling to extend the model to 3D. If I use spherical polar coordinates, $\phi$ is the angle between the Z-axis and the vector $\vec{OP}$, where P is the position of the planet.

As orbital inclination is measured from the positive x-axis to $\vec{OP}$, then $i=\frac{\pi}{2}-\phi$. From this diagram, $z=\frac{r}{tan(\phi)}$.

If I use the polar equation $r=\frac{a\left(1-e^2\right)}{1+ecos\ \theta}$, then $$z\left(t\right)=\frac{a\left(1-e^2\right)}{\left(1+ecost\right)tan\left(\frac{\pi}{2}-i\right)}=\frac{a\left(1-e^2\right)tan\left(i\right)}{\left(1+ecost\right)}$$

However, when I graph Venus and Earth's orbit on the same axis, I get this result in MATLAB. I don't believe there are any errors in my code, so it must be an error in my mathematics.

tl;dr Is there a better way to model ellipses that have a specific angle of inclination?

function [] = orbitEquations_3d(planetName, eccentricity, a, inclination)
% initialise (1-e^2) and sqrt(1-e^2)

    oneMinusE_sq = (1-eccentricity^2);
    sqrtOneMinusE_sq = sqrt(oneMinusE_sq);

    % convert angle of inclination from deg to radians
    inclination = inclination/180 * pi; 

% planet name
    disp(planetName);
    
% display equations
    disp("parametric");
    fprintf("x(t) = %.6f*cos(t) + %.6f\n", a, eccentricity*a);
    fprintf("y(t) = %.6f*sin(t)\n", a * sqrtOneMinusE_sq);
    fprintf("z(t) = %.6f*tan(%.6f)/(1+%.6fcos(t)) \n", a * oneMinusE_sq, inclination, eccentricity);
    fprintf("\n\n");
end