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




