Is it possible to create diff(function, variable) into a variable?

Hi! Thank you in advanced for the help.

I am trying build a script to solve a system of n pendulums. I am applying the same math as described here: Mathematical approach (I will further use some of the same notation and refer to this).

My goal is to end up with a set of expressions of ThetaDoubleDot_i, where i = 1,2,...,n, and these expressions should only depend on the variables Theta and ThetaDot.

I am using Euler-Langrange Differential Equation, which means that I need to be able to take the derivative of equation L (see URL) with respect to Theta_i, ThetaDot_i and time t.

MY APPROACH

I define the variables theta_i, thetaDot_i, theta_i(t), thetaDot_i(t) and t, where theta_i(t) and thetaDot_i(t) becomes functions of t. This enables me to first express L as expression of variables theta_i and thetaDot_i, and in that way easy use the diff-function in MatLab to solve its derivatives. Just before I take the derivate with respect to t (see eq. 11 and 16 in the attached URL) i substitute all theta_i and thetaDot_i with theta_i(t) and thetaDot_i(t), and this seems to work all fine.

PROBLEM

For instance one of the equations I end up with is ( w = Theta):

9.81*sin(w2(t)) + cos(w1(t) - w2(t))*diff(wDot1(t), t) + diff(wDot2(t), t) - sin(w1(t) - w2(t))*wDot1(t)*diff(w1(t), t) + sin(w1(t) - w2(t))*wDot1(t)*diff(w2(t), t) - sin(w1(t) - w2(t))*wDot1(t)*wDot2(t)

This is the same expression as eq. 18 in the URL, with alle lengths = 1, and masses = 1. But the problem now is:

I want to substitute back what I did earlier so that theta_i(t) becomes theta_i and so on. The reason for this is that then I will be able to handle the expression as a function of theta_i and thetaDot_i, instead of a function of t which it is while expressed by theta_i(t) and thetaDot_i(t). This is vital for the further calculations of my script. The problem is that i now have terms like diff(wDot1(t), t), which make the substitution wrong. Ideally I could convert diff(wDot1(t), t) into a variable ThetaDoubleDot_i (wDoubleDot1) and diff(w1(t), t) into thetaDot_i (wDot1). Is this possible?

My code so far:

    %%n pendulums 3.0
    % Sondre Nærland
    % 31. october 2018

    %{
    KOMMENTAR
    %}

    clear all
    close all
    clc

    %%Adjustable varables

    n = 2;                          % number of pendulums. Must be >=2.
    rod_length = ones(n,1);         % rod length is given the length of 1 for all rods. However adjustable.
    pointmasses = ones(n,1);        % masses of pointmasses. Given the value of 1 for all. However adjustable.
    anglesStart = ones(n,1);        % starting angle, given initial value of 1 for all. However adjustable.

    %%Calculated/fix variables

    g = 9.81;                        % gravity

    %%Math

    createdSymbols = symbolsFunc(n);  
    anglesTime = createdSymbols(:,1);
    anglesDotTime = createdSymbols(:,2);
    angles = createdSymbols(:,3);
    anglesDot = createdSymbols(:,4);
    syms t
    %{
    for i = 1:1:n
        sym(angles(i));
        sym(anglesDot(i));
    end
    %}
    xpos = sym('a', [n 1]);
    ypos = sym('b', [n 1]);
    xdot = sym('c', [n 1]);
    ydot = sym('d', [n 1]);
    xpos(1) = rod_length(1)*sin(angles(1));                            
    ypos(1) = -rod_length(1)*cos(angles(1));
    xdot(1) = anglesDot(1)*rod_length(1)*cos(angles(1));
    ydot(1) = anglesDot(1)*rod_length(1)*sin(angles(1));

    for i = 2:1:n
        xpos(i) = xpos(i-1) + rod_length(i)*sin(angles(i));
        ypos(i) = ypos(i-1) - rod_length(i)*cos(angles(i));
        xdot(i) = xdot(i-1) + anglesDot(i)*rod_length(i)*cos(angles(i));
        ydot(i) = ydot(i-1) + anglesDot(i)*rod_length(i)*sin(angles(i));
    end
    %%define L = Ek - Ep
    Ep = 0;
    Ek = 0;
    for i = 1:1:n
        Ep = Ep + pointmasses(i)*g*ypos(i);
        Ek = Ek + 0.5*pointmasses(i)*(xdot(i)^2 + ydot(i)^2);
    end
    L = Ek - Ep;
    % Ønsker d/dt(dL/d(generalized dotPosition)) - dL/d(generalized position) = 0
    listEq = sym('e', [n 1]);
    for i = 1:1:n
        d1 = diff(L, angles(i));    % dL/d(generalized position)
        d2 = diff(L, anglesDot(i)); % dL/d(generalized dotPosition)
        subD1 = subs(d1, [angles(:), anglesDot(:)], [anglesTime(:), anglesDotTime(:)]);
        subD2 = subs(d2, [angles(:), anglesDot(:)], [anglesTime(:), anglesDotTime(:)]);
        d3 = diff(subD2, t);           % d/dt(dL/d(generalized dotPosition))
        eq = d3 - subD1;
        listEq(i) = subs(eq, [anglesTime(:), anglesDotTime(:)], [angles(:), anglesDot(:)]);
    end

And the function that creates list of symbols

function list = symbolsFunc(n) 

list = sym('a',[n 4]);
for i = 1:1:n
    syms(sprintf('w%d', i))
    syms(sprintf('wDot%d', i));
    syms t
    w = symfun(sprintf('w%d(t)', i), t);
    wDot = symfun(sprintf('wDot%d(t)', i), t);
    wUT = symfun(sym(sprintf('w%d', i)), sym(sprintf('w%d', i)));
    wDotUT = symfun(sym(sprintf('wDot%d', i)), sym(sprintf('wDot%d', i)));
    list(i,1) = w;
    list(i,2) = wDot;
    list(i,3) = wUT;
    list(i,4) = wDotUT;
end