clear

L = 10;
EI = 100;
XY = [ ..
    1, 1;..
    2, 2;..
    3, -1..
];
springStiffness = 100;

function c = beamStiffness(beamLength)
    l = beamLength;
    l2 = l*l;
    l3 = l2*l;
    c = EI/l3 * [12,6*l,-12,6*l; 6*l,4*l2,-6*l,2*l2; ..
     -12,-6*l,12,-6*l; 6*l,2*l2,-6*l,4*l2];
endfunction

function u = bendingBySprings(xy, x)
    springCount = size(xy, 1);
    beamCount = springCount + 1;
    n = beamCount * 2;
    C = zeros(n,n);
    F = zeros(n,1);
    c = beamStiffness(xy(1,1));
    C(1:2, 1:2) = C(1:2, 1:2) + c(3:4, 3:4);
    for i=1:springCount
        xnext = L;
        if (i<springCount) then
            xnext = xy(i+1,1);
        end
        c = beamStiffness(xnext-xy(i,1));
        r = 2*i-1:2*i+2;
        C(r,r) = C(r,r) + c;
        r = 2*i-1;
        C(r,r) = C(r,r) + springStiffness;
        F(r) = F(r) + springStiffness*xy(i,2);
    end
    U = C\F;
    // TODO: Generate u
    u = x.*x;
endfunction

x = linspace(0, L, 500)';

u = bendingBySprings(XY, x);

plot(x, u);