Need to Smooth Plotted Curve

Good afternoon, folks. I've got an assignment I'm most of the way through but the final result is a little wonky looking. The two resulting plots should be relatively smooth curves but, while they have roughly the correct shape, the actual paths of the lines are weird in a way I'm not sure how to describe. I've attached both the problem statement and the data set necessary for the problem, and the code is:

function [t,x] = vib3
clear all
clc
global mp mr c k e w1 n w2 y_in an bn
% Define the parameters
mp = 300; % mass, kg, pump
mr = 240; % mass, kg, rotor
c = 2000; % damping constant, N-s/m
k = 200000; % stiffness, N/m
e = 0.02; % eccentricity distance, m
rpm = 3600;
w1 = (rpm * 2 * pi)/60; % forcing or input frequency, rad/s
wn = sqrt(k/mp);
zeta = c / (2 * mp * wn);
r = w1 / wn;
ratio = sqrt((1+(2*zeta*r))/((1-r^2)^2+(2*zeta*r)^2));

% Define the time interval for simulation
t_init = 0;
t_final = 5; % 0 to 5 sec.
t_inc = 0.01;
n = t_final/t_inc + 1; % 501 data points
t_span = linspace(t_init, t_final, n);
w2 = (2*pi)/t_final;

% Define the initial condition
x_init = 0.0; % initial position, m
xdot_init = 0.0; % initial velocity, m/s
x_0 = [ x_init; xdot_init ];

y = load('test_results.dat');
y_in = y(:,2) / 1000;
cn = fft(y_in)/length(y(:,1));
an = imag(cn);
bn = real(cn);



% Solve the system equations in differential form
[t,x] = ode45(@model, t_span, x_0);

% xdd = diff(x(:,2))/t_inc;
% ftrans=0;
% 
% for i = 1:length(xdd)
%     ftrans = ftrans + mp*xdd(i) - ratio*mr*e*w1^2*sin(w1*t);
% end

ftrans = k*x(:,1) + c*x(:,2);

subplot(1,2,1)
plot(t,x(:,1))
title('Displacement of TurboPump')
xlabel('time, t')
ylabel('Displacement, m')

subplot(1,2,2)
plot(t,ftrans/1000)
title('Force Transmitted to Base')
xlabel('time, t')
ylabel('Transmitted Force, kN')

end

function f = model(t,x)
global mp mr c k F w1 w2 n e y_in an bn
gral = trapz(y_in);
a0 = (w2 / pi)*gral;
y = a0 / 2;
ydot = 0;

for i = 1:n
    y = y + an(i)*cos(i*w2*t) + bn(i)*sin(i*w2*t);
    ydot = ydot - an(i)*i*w2*sin(i*w2*t) + bn(i)*i*w2*cos(i*w2*t);
end
F = k*y + c*ydot + mr*e*w1^2*sin(w1*t);
xdot1 = x(2);
xdot2 = ( F - k*x(1) - c*x(2) ) / mp;
f = [ xdot1; xdot2 ];
end

For context, we're meant to solve this using the fourier transform, which is why my y and ydot look like that.

I'm all but certain my error is in the setup of the for loop but I cannot figure out how I would change that. Also, if it isn't too much trouble, we're meant to be finding the ftrans value in the way I have commented out, but the plot resulting from that doesn't peter out the way the displacement curve does, and also gives values in the 10^8 N range, which seems high. Any advise on the matter would be greatly appreciated. Thank you!

[darova]

Can you show original/source formulas?

[Sarah Ten Eyck]

These are the notes we were given for the assignment, I hope they are what you are asking for.

[darova]

Are you sure i is a counter not a complex value?

for i = 1:n
    y = y + an(i)*cos(i*w2*t) + bn(i)*sin(i*w2*t);
%     y = y + an(i)*cos(1i*w2*t) + bn(i)*sin(1i*w2*t);
    ydot = ydot - an(i)*i*w2*sin(i*w2*t) + bn(i)*i*w2*cos(i*w2*t);
%     ydot = ydot - an(i)*1i*w2*sin(1i*w2*t) + bn(i)*1i*w2*cos(1i*w2*t);
end

[Sarah Ten Eyck]

I'm sure, yes. You can see, near the top of the notes, where I replace n with i.

[darova]

Theoretically: if t = 0 and n = inf, ydot = inf. BUt it's wrong. ISn't it?

ydot = ydot - an(i)*i*w2*sin(i*w2*t) + bn(i)*i*w2*cos(i*w2*t);

[Sarah Ten Eyck]

Sure, but we only take n to infinity in theory, otherwise it'd be functionally impossible to finish the equation, you'd never stop iterating. In reality you find a good enough value of n and just go with that.

[darova]

Run this script and see how y and ydot change. Look on the time interval [0 .. 0.2]

data = load('test_results.dat');
an = data(:,1);
bn = data(:,2);
w2 = 2*pi/5;
n = length(an);

for t = linspace(0,0.2,n);
    y = 0;
    dy = 0;
    for j = 1:n
        y = y + an(j)*cos(j*w2*t) + bn(j)*sin(j*w2*t);
        dy = dy - an(j)*j*w2*sin(j*w2*t) + bn(j)*j*w2*cos(j*w2*t);
    end
    subplot(211)
    hold on
    plot(t,y,'.r')
    subplot(212)
    hold on
    plot(t,dy,'.b')
end
hold off

[Sarah Ten Eyck]

Alright, I've done it, what am I meant to be getting from these plots?

[darova]

As you can see derivative ydot changes its sign many times. This might be the reason of why your main curve is not smooth