In order to generate an ellipse orientated along an axis that forms an angle φ with respect to the abscissa you have to apply the rotation of coordinates to the equation of the ellipse itself.
The equation of the ellipse is (as correctly implemented in you code):
xe=R•cos(θ) + xCenter
ye=R•sin(θ) + yCenter
the rotation is given by:
xr=xe•cos(φ)-ye•sin(φ)
yr=xe•sin(φ)+ye•cos(φ)
In order to have three ellipses 120° spaced, simply apply φ1=0°, φ2=120°, φ3=240°.
Since the axis rotation will rotate (of course) also the center, if xCenter=0 and yCenter=0 the centers af the three ellipses will be the same and equal to (0,0), so to have also a rotation of the center you will need to place it at a certain distance from the origin. To have the ellispses touching in a single point simply set to 0 the center corresponding to the minor axis of the ellipse and set the other at the same value of the radius of the major axis.
If, for instance: xRadius=8, yRadius=3, then xCenter=8, yCenter=0
and if xRadius=4, yRadius=7, then xCenter=0, yCenter=7
I wrote a code for Scilab. It is very similar to Matlab then with minor changes you can run it of that platform (for instance "isoview" stays for "axis square" and "%pi" stays for "pi").
xCenter = 0;
yCenter = 8;
xRadius = 2;
yRadius = 8;
theta = 0 : 0.01 : 2*%pi;
R=xRadius;
isoview(-20,20,-20,20);
x = xRadius * cos(theta)+xCenter;
y = yRadius * sin(theta)+yCenter;
phi=0*%pi/180;
xr = x*cos(phi)-y*sin(phi);
yr = x*sin(phi)+y*cos(phi);
plot(xr, yr, 'r');
phi=120*%pi/180;
xr = x*cos(phi)-y*sin(phi);
yr = x*sin(phi)+y*cos(phi);
plot(xr, yr, 'b');
phi=240*%pi/180;
xr = x*cos(phi)-y*sin(phi);
yr = x*sin(phi)+y*cos(phi);
plot(xr, yr, 'g');
If you need the ellipses to be intersecated, simply decrease the value of yCenter.