solving beam deflection curve

I had a question about plotting Beam deflection curve

The drawing was almost correct, but i had a percentage of error when measuring the deflection at 1500 and 3000 compared to the results I got from skyciv.

I do not understand why this break occurred in the third period, and despite my application of the boundary condition equation at 2000 that the equation value equal zero, but in the loop and on the curve the value is not zero at the support third support

clc,clear
E = 250e3; % Young's modulus
I = 2e7; % Moment of inertia
L = 1000; % Length of the beam
F = 5000; % Force at A
W1 = 200; % Distributed load 1
W2 = 20; % Distributed load 2
W3 = 0.1; % Distributed load 3
m = W3*(1000^2)/2

m = 50000

syms R1 R2 R3 c
eqf = R1 +R2+R3 - F -W1*L/3-W3*L-W2*L ==0

eqf = 

eqm = R2*L+2*R3*L-F*L/3-W1*L/3*0.5*L-W3*L*1.5*L+W3*L^2/2-2.5*W2*L^2==0

eqm = 

cdef = (R1*L^3)/6-(4/81)*F*L^3-W1*L^4/144+c*L == 0

cdef = 

def2 = R1*4*L^3/3+R2*L^3/6-(125/162)*F*L^3-3/16*W1*L^4-W3*L^4/48-W3*L^4/16+2*L*c==0

def2 = 

sol = solve(eqf,eqm,cdef,def2)

sol = struct with fields:
    R1: 2183306111051851699891/67108864000000000
    R2: 1119568478814815206359/33554432000000000
    R3: 3471827368770369759157/134217728000000000
     c: -1524645038459259098803/402653184000

R1 = double(sol.R1)

R1 = 3.2534e+04

R2 = double(sol.R2)

R2 = 3.3366e+04

R3 = double (sol.R3)

R3 = 2.5867e+04

c = double(sol.c)

c = -3.7865e+09

x = linspace(0, 3000, 300001);
y = zeros(size(x));
% Compute y values
for n = 1:length(x)
    if x(n) < (1000/3)
        y(n) = (1/(E*I)) * ((R1*x(n)^3)/6 + c*x(n));
    elseif x(n) >= (1000/3) && x(n) < (2000/3)
        y(n) = (1/(E*I)) * ((R1*x(n)^3)/6 - F*(x(n)-(L/3))^3/6 - W1*(x(n)-(L/3))^4/24 + c*x(n));
    elseif x(n) >= (2000/3) && x(n) < 1000
        y(n) = (1/(E*I)) * ((R1*x(n)^3)/6 - F*(x(n)-(L/3))^3/6 - W1*L/3*(x(n)-(L/2))^3/6 + c*x(n));
    elseif x(n) >= 1000 && x(n) < 1500
        y(n) = (1/(E*I)) * ((R1*x(n)^3)/6 - F*(x(n)-(L/3))^3/6 - W1*L/3*(x(n)-(L/2))^3/6 + R2*(x(n)-L)^3/6 + c*x(n));
    elseif x(n) >= 1500 && x(n) < 2000
        y(n) = (1/(E*I)) * ((R1*x(n)^3)/6 - F*(x(n)-(L/3))^3/6 - W1*L/3*(x(n)-(L/2))^3/6 + R2*(x(n)-L)^3/6 - (W3*((x(n)-3/2*L)^3)/6 + m*((x(n)-(3/2)*L)^2)/2) + c*x(n));
    elseif x(n) >= 2000 && x(n) <= 3000
        y(n) = (1/(E*I)) * ((R1*x(n)^3)/6 - F*(x(n)-(L/3))^3/6 - W1*L/3*(x(n)-(L/2))^3/6 + R2*(x(n)-L)^3/6 - (W3*((x(n)-3/2*L)^3)/6 + m*((x(n)-(3/2)*L)^2)/2) + R3*((x(n)-2*L)^3)/6 - W2*((x(n)-2*L)^4)/24 + c*x(n));
    end
end
plot(x,y)
grid on

YC  = (1/(E*I)) * ((R1*1500^3)/6 - F*(1500-(L/3))^3/6 - W1*L/3*(1500-(L/2))^3/6 + R2*(1500-L)^3/6 + c*1500)

YC = 0.1762

YE =  (1/(E*I)) * ((R1*3000^3)/6 - F*(3000-(L/3))^3/6 - W1*L/3*(3000-(L/2))^3/6 + R2*(3000-L)^3/6 - (W3*((3000-3/2*L)^3)/6 + m*((3000-(3/2)*L)^2)/2) + R3*((3000-2*L)^3)/6 - W2*((3000-2*L)^4)/24 + c*3000)

YE = -1.2924

Y_C = (1/(E*I))*((R1*3*L^2)/2 -F*((3*L-(L/3))^2)/2 - W1*L/3*(3*L-(L/2)^2)/2+(R2*(3*L-L)^2)/2+c)

Y_C = 0.0204

Y_E = (1/(E*I))*((R1*3*L^2)/2 -F*((3*L-(L/3))^2)/2 - W1*L/3*(3*L-(L/2)^2)/2+(R2*(3*L-L)^2)/2-(W3*((3*L-3/2*L)^2)/2+ m*(3*L-(3/2)*L))+(R3*((3*L-2*L)^2)/2-W2*(((3*L-2*L)^3)/6))+c)

Y_E = 0.0223

Answer 1

[VBBV]

if you use linspace function, it divides the length between initial and final points into non-equispaced values, however if you specify vector of values between 0 and 3000, then it does consider each value which includes points 1500 and 3000 of the beam.

x = 0:1:3000;
Ix = x(x==1500)  % find the index of location 1500 

Ix = 1500

x = linspace(0,3000,30000)

x = 1×30000

0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000 1.1000 1.2000 1.3000 1.4000 1.5001 1.6001 1.7001 1.8001 1.9001 2.0001 2.1001 2.2001 2.3001 2.4001 2.5001 2.6001 2.7001 2.8001 2.9001

Ix = x(x==1500)  % no such point exists in the vector when you use linspace function

Ix =

  1×0 empty double row vector