I ended up doing a mixture of solutions. Thanks @Marco, @Cris Luengo from stackoverflow and @Karim from matlab's forum.
- r is the camera position relative to the origin like everything else
- dist is half the distance of one of the sides of the cube displayed. I decided to use the orbital height multiplied by tan(68º), based on human horizontal FOV
- rSun is the position vector of the sun
- rMoon (drawing this one too) is the position vector of the moon
First of all I have an infinite light source with the position of the sun.
To draw a far away surface object like the moon I place it with a special rMoon2 that locates it in the edge of the plot and multiply its size by a proportional number:
rMoon2 = rMoon/norm(rMoon)*dist
size = size*(dist-norm(r))/(norm(rMoon)-norm(r))
If the object is extremely far away like the sun you can drop the denominator's r.
If the object is itself a light source (the sun), you need another light source (but local) between the observer located in r, and the sun. So additionally to earlier steps you need:
localLightPos = rSun2*(1-2*tand(34/60)*(dist-norm(r))/dist))
Where 2*tand(34/60)*(dist-norm(r)/dist) is twice the radius of the apparent sun relative to the distance between sun and observer.
Here is the entire code for the plot in case you want it.
r = Y(j,1:3);
dist = tand(68)*(norm(r)-Rearth);
figure('Color','k','Position');
% Celestial bodies
p3Dopts.Units = 'km';
p3Dopts.RotAngle = angleEarth(j);
planet3D('Earth', p3Dopts);
hold on
p3Dopts = rmfield(p3Dopts, 'RotAngle');
p3Dopts.Position = rSun(j,:)/norm(rSun(j,:))*dist';
relDist = 2.1*(dist-norm(r)); % 2.1 because the sun appears bigger than it is
p3Dopts.Size = relDist/norm(rSun(j,:));
planet3D('Sun', p3Dopts);
%Sunlight, apparent Sun is ~28-34 arc minutes near Earth
light("Style","infinite","Position",p3Dopts.Position)
light("Style","local","Position",p3Dopts.Position*(1-2*tand(34/60)*relDist/dist));
p3Dopts.Position = rMoon(j,:)/norm(rMoon(j,:))*dist';
p3Dopts.Size = relDist/norm(rMoon(j,:)-r);
planet3D('Moon', p3Dopts);
p3Dopts = rmfield(p3Dopts, 'Position');
p3Dopts = rmfield(p3Dopts, 'Size');
% Ground track
scatter3( trackT(1:j,1), trackT(1:j,2), trackT(1:j,3), 12, T(1:j), 'filled')
%xlabel('x [km]'); ylabel('y [km]'); zlabel('z [km]');
title('Earth focused animation');
axis equal;
grid on;
ax = gca; ax.Color = 'k'; ax.GridColor = 'k';
ax.GridAlpha = 0; ax.XColor = 'k'; ax.YColor = 'k';
ax.ZColor = 'k';
ax.CameraPosition = r; % Set the camera position
ax.XLim = [-dist, dist]; % Set the x-axis range
ax.YLim = [-dist, dist]; % Set the y-axis range
ax.ZLim = [-dist, dist]; % Set the z-axis range