Point by Point Method

pbp.m

clc;
clearvars;

t = 0;
tf = 0;
tfinal = 0.5;
tc = 0.125;
tstep = 0.005;
M = 2.52/(180*50);
i = 2;
delta = 21.64*pi/180;
ddelta = 0;
time(1) = 0;
ang(1) = 21.64;
Pm = 0.9;
Pmaxbf = 2.44;
Pmaxdf = 0.88;
Pmaxaf = 2.00;

while t < tfinal
   if t==tf
       Paminus = 0.9 - Pmaxbf*sin(delta);
       Paplus = 0.9 - Pmaxdf*sin(delta);
       Paav = mean([Paminus Paplus]);
       Pa = Paav;
   end
   if t==tc
       Paminus = 0.9 - Pmaxdf*sin(delta);
       Paplus = 0.9 - Pmaxaf*sin(delta);
       Paav = mean([Paminus Paplus]);
       Pa = Paav;
   end
   if t>tf && t<tc
       Pa = Pm - Pmaxdf*sin(delta);
   end
   if t>tc
       Pa = Pm - Pmaxdf*sin(delta);
   end
   
   ddelta = ddelta + (tstep^2*Pa/M);
   delta = (delta*180/pi + ddelta)*pi/180;
   deltadeg = delta*180/pi;
   time(i) = t;
   ang(i) = deltadeg;
   t = t+tstep;
   i = i+1;
end
figure(1);
axis([0 0.6 0 160]);
xlabel('Time');
ylabel('Angle');
plot(time,ang, 'k-');

Some change is needed, for the swing angle curve; plot of given code is not a swing eq. plot.

clc;
clearvars;

t = 0;
tf = 0;
tfinal = 0.5;
tc = 0.125;
tstep = 0.005;
M = 2.52/(180*50);
i = 2;
delta = 21.64*pi/180;
ddelta = 0;
time(1) = 0;
ang(1) = 21.64;
Pm = 0.9;
Pmaxbf = 2.44;
Pmaxdf = 0.88;
Pmaxaf = 2.00;

while t < tfinal
   if t==tf
       Paminus = 0.9 - Pmaxbf.*sin(delta);
       Paplus = 0.9 - Pmaxdf.*sin(delta);
       Paav = mean([Paminus Paplus]);
       Pa = Paav;
   end
   if t==tc
       Paminus = 0.9 - Pmaxdf*sin(delta);
       Paplus = 0.9 - Pmaxaf*sin(delta);
       Paav = mean([Paminus Paplus]);
       Pa = Paav;
   end
   if t>tf && t<tc
       Pa = Pm - Pmaxdf*sin(delta);
   end
   if t>tc
       Pa = Pm - Pmaxaf*sin(delta);
   end
   
   ddelta = ddelta + (tstep^2*Pa/M);
   delta = (delta*180/pi + ddelta)*pi/180;
   deltadeg = delta*180/pi;
   time(i) = t;
   ang(i) = deltadeg;
   t = t+tstep;
   i = i+1;
end

figure('Name', 'Swing Eq. when Fault Sustain')
axis([0 0.6 0 160]);
p = plot(time,ang, 'color', 'r', 'LineWidth', 2);
xlabel('Time');
ylabel('Angle');
title('Swing Eq. when Fault Sustain'); grid on;
str = {['Maximum Angle is ', num2str(max(ang)), ' at time of ', num2str(time(find(ang == max(ang))))]};
annotation('textbox',[0.25 0.1 0.3 0.2],'String',str,'FitBoxToText','on');
datatip(p,time(find(ang == max(ang))) , max(ang))