% Use Symbolic Toolbox to perform integration
syms P L EI x k positive;

M = sym('P*x');             % M equation
EI = sym('EI*(1+k*x/L)');   % EI equation
Mohr2 = M/EI*x;             % 1st moment of M/EI diagram
def = int(Mohr2,x,0,L);     % definite integral
pretty(def);                % display result

limit(def,k,0);             % Prove limit as k->0 is prismatic result

% Plot change in deflection by varying k
clear all;
P = 10; L = 4; EI = 10e3; k = 0.01:0.01:2;
d = 1/2*P*L^3*(-2.*k+2*log(1+k)+k.^2)/EI./k.^3;
plot(k,d*1e3); xlabel('Stiffness Increase at A (k)'); 
ylabel('Deflection (mm)')

d1 = P*L^3/(3*EI);          % prismatic result