The following is a very, very simple program modelling a variable star, using octave. There's very little physics in it. t=0 r=1 ampl=0.01 dt=5 # minutes period= 30.0 for i=2:100 t = [t, t(i-1)+dt]; r = [r, r(1) + ampl * sin( pi*( ((t(i)-t(1)) / period ) ) ) ] ; endfor plot (t,r)