You've only got four points there, when you need four orbits, right?
I don't know enough about equations for elliptical orbits to know where your theta are supposed to fit in, but I know that you'll have to specify points all the way around the orbit:
p = [.387 .723 1 1.524];
E = [.2056 .0068 .0167 .0934];
% theta = deg2rad([7.005 3.3947 0 1.851]);
theta = linspace(0,2*pi,100).';
radius = p./(1 + (E.*cos(theta)));
polarplot(theta, radius,'LineWidth',2);
That gives you some orbits; maybe you can work with it to get what you're going for.
