Verify the solution of beam loading given a sinusodial function

%% I need help implementing this specific part (part 2 and later part 3): Where we have a stiffness matrix of a beam and I need to plot the load of the given function. I have produced an graph but it looks odd, so im not sure if i implemented the code correctly (Part 2)

.

=============================================================================

%% main script
%% Constants
L=2;
w=.3;
d=.03;
E=13e3;
rho=480;
g=9.81;
I=  (w*d^3)/(12);
%% Part 1 create a NxN matrix of K using stiff_m    atrix function
N = 250;%input('please enter a number N to create the stiffness matrix: ');
h=L/N;
k_coeff = (E*I)/(h^4);
k = k_coeff* stiff_mat(N);
%% Part 2: Sinusodial Loading
x=linspace (0,N,N);
L=2;
fn = sin_func(x);
% calculate y(displacement) =k/fn
y=fn.'\k;
figure
hold on
plot([0 N],[0 0],'k-'); %plots flat line (y=0)
plot(x,y,'k-'); %plots curve of y
xlim = ([0,N]);
title('graph of loading')

%% Part 3: Validate the Method
y_x= ((-1*rho*g*w*d)/(24*E*I)) .*(x.^2).*6*L^2 - 4*L.*x + x.^2;
y_L = ((-1*rho*g*w*d*L^4)/(8*E*I));

%==============================================================
%% Define stiffness matrix for a given NxN dimension
function [stiff_matrix] = stiff_mat(N)
stiff_matrix = zeros(N);
stiff_matrix = stiff_matrix+6*eye(N);
minus_fours_above = diag(-4*ones(N-1,1),1);
minus_fours_below = diag(-4*ones(N-1,1),-1);
ones_below = diag(ones(N-2,1),-2);
ones_above = diag(ones(N-2,1),2);
stiff_matrix = stiff_matrix +minus_fours_above + minus_fours_below +ones_above+ones_below;
stiff_matrix(1,1) = stiff_matrix(1,1)+1;
stiff_matrix(1,3) = stiff_matrix(1,3)+2;
stiff_matrix(N-1,N-1) = stiff_matrix(N-1,N-1)-1;
stiff_matrix(N,N) = stiff_matrix(N,N)-4;
stiff_matrix(N-1,N) = stiff_matrix(N-1,N)+2;
stiff_matrix(N,N-2) = stiff_matrix(N,N-2)+1 ;
end
%============================================================
%% sinusodial function
function [y] = sin_func(x)
L=2;
y = 100.*sin((4*pi.*x)/(3*L));
end

[Dyuman Joshi]

The expression for y_x is incorrect -

y_x= ((-1*rho*g*w*d)/(24*E*I)) .*(x.^2).*6*L^2 - 4*L.*x + x.^2;

it should be

y_x= ((-1*rho*g*w*d)/(24*E*I)) .*(x.^2).*(6*L^2 - 4*L.*x + x.^2);

And you can use the sin_func directly (L is already defined in the code)

fn = 100.*sin((4*pi.*x)/(3*L));

Rest looks good.