I don't think the results are going to be very interesting but you could do what you are asking as follows:
ro=1.49e11; %Sun orbital radius (1 AU)
G=6.67408e-11;
MS=2e+30;
%synodic period
perp = 365; % planet period (days)
perd = 1896.59; % red dwarf period (days)
np = 2*pi / perp; % mean motion of planet (rad/day)
nd = 2*pi / perd; % mean motion of red dwarf (rad/day)
dT = 2*pi / (np - nd); % This is the average time between stellar eclipses.
v=sqrt((G*MS)/ro);
T = dT:dT:dT*ceil(100*perp/dT); % up to next multiple of dT after 100 years (multiples of perp)
theta=((v*T)*86400)/ro % Angle in radians at each multiple of dT
This is just going to add a constant increment to theta every dT. If you plot(T, theta) you will get a boring straight line. You could, of course, apply a mod by 2*pi but this won't be much more useful (IMHO).
Anyway, I hope this answers your question. If not, please clarify.