ClearAll["Global`*"]
rstar=2(6.957*10^8)​​pts=50​​dx=2*rstar/pts;​​dy=2*rstar/pts;​​S=Table[(1-0.5(1-(3/2)Sqrt[1-(Sqrt[(x)^2+(y)^2]/rstar)^2])),{x,-rstar,rstar,dx},{y,-rstar,rstar,dy}];​​MatrixForm[S]​​ListDensityPlot[S,DataRange{{-rstar,rstar},{-rstar,rstar}},PlotLegendsAutomatic,ColorFunction"SolarColors"]
1.3914×
9
10
50