clear all
syms a x n_terms v PI EE II L q
figure
for n_terms=3:6 % Try solutions between 3 and 6 terms
    clear eq
    a = sym('A',[n_terms,1]); % Create a vector containing the unknown coefficients
    v = 0; % Define the approximation
    for i=1:n_terms
        v=v+a(i)*x^(i-1);  %Adds each term to the approximation one at a time
    end
    n_constraints = 2;  % Enforce the constraints - v and dv/dx=0 at x=0
    eq(1) = subs(v,x,0)==0;
    eq(2) = subs(diff(v,x),x,0)==0;
    asol = struct2cell(solve(eq,a(1:n_constraints))); % Eliminate two of the unknowns 
    for i=1:n_constraints
        v = subs(v,a(i),asol{i});
    end
    PI = int(EE*II/2*diff(v,x,2)^2 - q*v,x,[0,L]); % Potential energy
    for i=n_constraints+1:n_terms
        eq(i)=diff(PI,a(i))==0;  % Differentiate wrt all the remaining coefficients
    end
    if (n_terms>n_constraints+1) % Solve for the remaining unknowns
       asol = struct2cell(solve(eq(n_constraints+1:n_terms),a(n_constraints+1:n_terms)));
       for i=n_constraints+1:n_terms
         v = subs(v,a(i),asol{i-n_constraints});
       end
    else
        asol = solve(eq(n_constraints+1),a(n_constraints+1));
        v = subs(v,a(n_constraints+1),asol);
    end
    v % Display the results and plot it
    string = [num2str(n_terms) ' Terms'];
    fplot(subs(v,[EE,II,L,q],[1,1,1,1]),[0,1],'DisplayName',string)
    hold on
end
axes1 = gca;
ylabel({'vEI/qL^4'});
xlabel({'x/L'});
box(axes1,'on');
axis(axes1,'tight');
% Set the remaining axes properties
set(axes1,'FontSize',18,'XGrid','on','YGrid','on');
legend1 = legend(axes1,'show');
set(legend1,...
    'Position',[0.245245334727259 0.532128183414017 0.233928567171097 0.270883047210857]);